Mathematical Model of the Decomposition of Unstable Gas Hydrate Accumulations in the Cryolithozone

Abstract

We present a generalization of the mathematical model of gas discharge from frozen rocks containing gas-saturated ice and gas hydrates in a metastable state (due to the self-preservation effect) caused by the drop in external stress associated with various geodynamic factors. These factors can be attributed, for example, to a decrease in hydrostatic pressure on a gas-bearing formation due to glacier melting, causing an isostatic rise, or to the formation of linear depressions in the bottom topography on the shelf due to iceberg ploughing. A change in external pressure can also be associated with seismic and tectonic deformation waves propagating in the lithosphere as a result of ongoing strong earthquakes. Starting from the existing hydrate destruction model, operating at the scale of individual granules, we consider a low-permeable hydrate and ice-saturated horizontal reservoir. Generalization is associated with the introduction of a finite threshold for the external pressure drop, which causes the destruction of the gas hydrate and gas-saturated microcavities of supramolecular size. This makes it possible to take into account the effect of anomalously high pressures occurring in the released gas as a result of partial hydrate dissociation. Numerical and approximate analytical solutions to the problem were found in the self-similar formulation. A parametric study of the solution was carried out, and regularities of the hydrate decomposition process were revealed.

1. Introduction

According to present-day concepts, gas hydrates, which are crystalline compounds of gas and water, are frequently found in the sedimentary rocks in permafrost areas both on land and in deep seas and oceans, where relatively high pressures and low temperatures occur, which are necessary for the hydrate formation and stable existence. At the same time, gas hydrates often occur at shallow depths under conditions of low pressures and negative temperatures, being in a metastable state due to the so-called self-preservation effect [1,2,3]. Existing estimates suggest that the global reserves of gas hydrates count for about 1000 billion metric tons of carbon [4], and the decomposition of such an amount can affect the carbon cycle and climate on a global scale [5,6,7,8,9]. The idea of possible gas hydrate dissociation causing methane emissions contributing to the greenhouse effect has been used to explain a number of well-known phenomena, such as the Paleocene-Eocene Thermal Maximum [5], or a rapid post-glacial increase in atmospheric methane concentration [10]. A general hypothesis is also discussed about the possible impact of methane emission bursts on global warming, with an increase in permafrost degradation and gas hydrate dissociation [9,11,12,13,14,15,16,17].

Long-term marine studies of the Russian Arctic Shelf have shown that a significant methane release occurs from the bottom of shallow water areas in the Eastern Arctic Seas [15,16], where the amount of released gas is comparable or exceeds the estimates accepted for the entire World Ocean [8,18]. Among numerous onshore gas manifestations of natural and man-made origin in the Arctic zone, a special place is occupied by large craters resulting from powerful gas emissions [19,20,21,22]. The formation of such craters is a potential hazard to human life, especially near the oil and gas field facilities. In addition, large methane seeps in the seas of the Eastern Arctic, which can pose a potential danger in the Northern Sea Route region. Mathematical modeling of heat and mass transfer in gases and fluids in the sedimentary formations containing frozen rocks and gas hydrates is a difficult problem, primarily due to the limited knowledge of the physical–chemical and mechanical properties of a multicomponent medium characterized by the gas hydrates’ dissociation and fluid filtration. It also involves phase transitions in the ice–water–gas hydrate system occurring inside the fracture pore space of the heterogeneous mineral matrix of host sedimentary rock. Usually, the gas hydrate decomposition is associated with an increase in temperature to a critical level, causing the loss of the hydrate stability at a given ambient pressure, corresponding to a certain depth of the hydrate layer. Such an increase in the temperature of the hydrate layer can occur due to either endogenous processes, e.g., when the heated fluid migrates along crustal faults from deeper parts of the lithosphere [23], or exogenous factors, such as warm near-bottom currents on the shelf or water transgressions on the cold land surface in the Arctic [12,13,15,18,24]

Atmospheric pressure fluctuations play an important role in determining the timing and quantity of these emissions [25], while tidal pressure acts as an important control mechanism on bubble release in shallow estuarine and marine environments [26]. Low air pressure events were associated with a storm system which appeared to induce ebullition from the northern lakes, whereas high pressure inhibited ebullition. Sometimes, the CH₄ ebullition increased up to two orders from “background values” and achieved values of tens of mM CH₄ m⁻² per day, when air pressure dropped to 740 Torr [27]. Thus, during the period of air pressure drop, the net CH₄ flux to the atmosphere was increased significantly. Such a correlation was obtained at Mirror Lake, New Hampshire [28]. A similar phenomenon had been known to colliery ventilation engineers in the United Kingdom for more than 250 years [29]. For example, the explosion at Haswell colliery, Durham, in 1844, was investigated by Michael Faraday and Charles Lyell. Their report suggested that a fall in the barometer resulted in the emission of methane into the airstream of the mine.

However, the loss of gas hydrate stability can also be attributed to the decrease in external pressure caused by geodynamic factors. Such factors include, for example, a decrease in the hydrostatic water pressure on the shelf as a result of its shallowing caused by the isostatic rise of the crustal surface due to glacier melting, as it was observed in the Svalbard archipelago in the Western Arctic [30]. Another exogenous geodynamic factor may be related to the widespread linear depressions in the bottom topography on the shelf, which are formed as a result of its ploughing by underwater keels of moving icebergs, known as gouging [31]. One of the clearest examples of massive methane discharge and the pockmark formation in the Barents Sea is the study [32], which substantiates the hydrate destabilization as a result of the sheet glacier destruction during the Holocene. Geodynamic factors can also play an important triggering role in the process of disturbing the metastable equilibrium of relict gas hydrates that experienced partial dissociation, which stopped as a result of the formation of ice films sealing the free gas inside gas hydrate microparticles, known as hydrate self-preservation [1,2,3,33]. Under self-preservation conditions, the metastable gas hydrates can break down when additional shear stresses occur, which may destroy thin ice films and cause free gas release, filtration and, ultimately, emission. Geodynamics-related causes for the occurrence of additional shear stresses in the lithosphere and sedimentary formations can be associated with seismic and tectonic deformation waves generated by earthquakes. Tectonic deformation waves produced by the largest subduction earthquakes and propagating in the elastic lithosphere–viscous asthenosphere system at a speed of several tens of kilometers per year are of particular importance [34]. Affecting vast areas (measuring over thousand kilometers), these rather slow waves bring additional stresses into the region where metastable gas hydrates are hosted, creating favorable conditions for the ice cover destruction and the release of methane trapped in micropores. This seismogenic triggering mechanism was proposed as a hypothesis linking methane emission bursts in the Arctic and the onset of recent climate warming phases in the 1920s and 1980s as a result of anomalously high seismic activity in the Aleutian subduction zone observed approximately 20 years prior to each phase at the beginning and in the middle of the 20th century [17,35].

It should be noted that in addition to the above, gas hydrates are involved in a variety of applications, which include production, transportation and storage of the gas hydrates; CO₂ absorption; water desalination technology; gas hydrate combustion technologies, etc. (see [36] and references therein). The model developed in this study can also be useful for research in those areas.

Ref. [37] studied the mechanism of the micropore-trapped gas filtration associated with the destruction of the pore-surrounding material (ice) due to the drop in external pressure. Their study analyzed the medium with agglomerate inclusions consisting of gas hydrate, ice and free gas. Although the porosity of such inclusions may be relatively high, they have extremely low permeability of the order of nanodarcy (10⁻²¹ m²) due to the nanosize of pores and channels. The model developed in [37] assumes that the effective permeability of the hydrate-bearing rock in which the destruction occurs is proportional to the pressure gradient to the power $m\ge 0$.

In this study, we develop the gas discharge model initially proposed in [37]. In particular, we avoid the somewhat artificial assumption made in the original model, stating that any small drop in external pressure causes the gas discharge from inclusions due to very small pore size. We believe that, in the general case, a critical pressure drop value exists, which may not be small, and depends on temperature, texture, and other factors, and which triggers the hydrate depreservation. This generalization is significant, since it allows us to take into account the effect of anomalously high pressure in the gas released from hydrates. In addition, the model discussed in this paper deals with a horizontal reservoir of abnormally low permeability, rather than individual inclusions. Thus, the original model [37] is the limiting case of our model, when the critical pressure differs very little from the unperturbed pressure in the reservoir. For gas hydrates, the unperturbed pressure is understood as the latent (potential) pressure of the gas, which exists during the hydrate decomposition in the absence of filtration. Under the conditions in question, the rock temperature usually changes insignificantly; therefore, we apply isothermal approximation. In the future study, it is planned to generalize the results to the nonisothermal case.

The most important factor controlling the hydrocarbon filtration is the rock permeability. In this regard, in conclusion of this section, we review a number of studies focusing on this parameter, primarily related to low-permeability media.

The study [38] discovered that low-permeability hydrate reservoirs have great potential for gas production claims, which was confirmed by the successful gas production tests in the South China Sea in 2017. To improve the energy efficiency of future commercial hydrate production, predict gas/water production rates and evaluate the multiphysics response of low-permeability reservoirs to gas production caused by depressurization, the study employs a fully coupled heat transfer and hydrodynamical simulation involving the Mohr–Coulomb geomechanical model. The above model is used to systematically study hydrate reservoir response in terms of pore pressure, temperature and hydrate saturation, and evaluate the mechanical behavior of a hydrate reservoir during one year of gas production using a horizontal wellbore. The simulation results show that the gas yield increases at the beginning of production, then quickly decreases after the cessation of further decrease in well pressure, and finally remains at the level of 100 m³/day. Due to the low permeability, the decreased pressure region extends up to 100 m, which is not very far from the wellbore. One year of gas production by a horizontal wellbore does not lead to a serious formation displacement. However, stress relaxation due to hydrate dissociation can cause a stress redistribution, leading to a displacement of about 0.35 m around the wellbore, which can potentially cause the well damage, affecting the safety and efficiency of gas production.

Paper [39] concludes that the macroscopic parameters of oil migration processes are controlled by the microscopic properties of the pore space in low-permeability media. To analyze the relationship between the pore scale parameters and the process of oil migration in low-permeability reservoirs, a mathematical model was developed in the above study that combines the capillary bundle model and the fractal theory. The accuracy of the proposed model is verified by comparison with three well-designed oil loading experiments using core samples. Based on the validated model, the influence of four factors (fractal dimensions of the pore space and tortuosity, rock surface wettability and pore water viscosity) on oil migration processes in low-permeability media was investigated. Oil saturation and effective permeability were used as output parameters indicating the oil migration. The calculation results show that oil saturation decreases with the increase in the pore space and tortuosity fractal dimensions, rock hydrophilicity and the brine viscosity. On the other hand, the effective permeability decreases with an increase in the fractal dimension of the tortuosity and pore water viscosity but increases with an increase in the pore space fractal dimension. In addition, the dependence of effective permeability on wettability was found to be relatively weak.

Work [40] applied nuclear magnetic resonance (NMR) for indirect detection and analysis of methane hydrate in subsurface formations, which is potentially useful for quantifying pore size in hydrate reservoir and for evaluating the hydraulic permeability of hydrate-bearing formations. In this experiment, hydrates were artificially formed on the seafloor by injecting methane into pipes with seawater-containing sediments; a few weeks later, the samples were examined remotely with an NMR instrument. A decrease in the NMR signal was associated with the hydrate formation, and the hydrate content determined by NMR allowed for predicting the volume of gaseous methane released during subsequent dissociation.

Studies [41,42,43] have shown that the water retention curve and relative permeability are critical for predicting gas and water production from hydrate deposits. However, the values of the key parameters characterizing the gas and water flow during the hydrate dissociation are hard to determine due to experimental problems. The study uses combined microfocus X-ray computed tomography (CT) and pore network modeling techniques to find the correct values for these key parameters, such as gas inlet pressure and residual water saturation. Using the micron-resolution CT images, a pore space geometry and saturation were recovered, and further used to simulate gas flow, hydrate dissociation and fluid permeability. The calculations show that higher hydrate saturation results in higher gas inlet pressure, higher residual water saturation and a steeper water retention curve. An increase in hydrate saturation reduces gas permeability but has little effect on water permeability in sediments with evenly distributed hydrate, while hydrate morphology has a more significant effect on relative permeability. Heterogeneously distributed hydrate accumulations, as a rule, have lower residual water saturation and higher gas and water permeability. In this sense, the Brooks–Corey model, using two appropriate parameters separately for gas and water permeability, correctly reflects the effect of hydrate saturation and morphology on gas and water flows in hydrate-bearing deposits.

Study [44] combined in situ microfocus X-ray computed tomography (CT) and pore network model simulation to obtain information on the dynamic permeability of hydrate dissociation, acknowledging that permeability is the most important factor controlling the formation of water and gas in hydrate sediments. The paper [45] emphasizes the importance of the pore hydrate phase equilibrium studies, including the determination of the equilibrium pore water content (non-clathrated water). Non-clathrated water is similar to unfrozen water in permafrost and has a significant impact on the hydrate-bearing reservoir properties, including permeability.

2. Materials and Methods

2.1. The Mechanism of the Destruction of Metastable Gas Hydrates

As is known, the decomposition of the hydrate might be caused by either an increase in temperature or a decrease in external pressure, and in this study, we focus on the second cause. Following modern idealized concepts, the gas hydrate in the boundary layer can be approximately viewed as a set of fine spherical granules (Figure 1a) [46]. Then, the mechanism of the self-preserved gas hydrate destruction can be analyzed on a scale of a single multilayer particle. We assume the simplest model of a spherical particle consisting of three spherical layers with different thicknesses (Figure 1b) [46,47]. The outer layer (1) is an ice crust formed during the gas discharge following the hydrate stability violation. Beneath the ice layer, there is a free gas layer underlain by the undecomposed hydrate spherical core. The pressure in the free gas is no lower than the pressure of the phase equilibrium of the hydrate, and therefore, this prevents the decomposition of the gas hydrate occupying the core region (3). Other existing models are essentially similar to the one described above. For example, the model proposed in [46] (Figure 1c) has no isolated gas layer, instead it is assumed that the ice crust with a thickness ${\delta }_2$ contains gas-filled pores. When the phase transition front (Ri) moves, the gas leaves the granule through the pores, formed as a result of a phase transition since the hydrate’s density is lower than that of the ice, and the empty microcavities are filled with gas. In the framework of such models, the problem of gas hydrate decomposition in the boundary layer is reduced to the problem of the ice crust breakdown in individual particles, and evaluation of the effect of temperature and pressure drop, etc., on the strength of the ice layer.

In [47], authors evaluate strength and breakdown conditions for a gas hydrate granule assuming the model qualitatively similar to that shown in Figure 1b. Their estimates are based on experimental data on the tensile strength of ice depending on the grain size obtained in [48] for temperature of −10°C and strain rate of 10⁻⁶ s⁻¹. Shimada et al. assume the granule geometry in which the gas layer (2) is thin, so that its thickness can be neglected, the ice layer (1) is 0.5 mm thick and the region (3) occupied by the undecomposed gas hydrate has a radius of 2 mm. The temperature is −10°C, the gas pressure is equal to the phase equilibrium pressure and the gas pressure outside the granule is equal to the atmospheric pressure. Estimates have shown that for an ice crust with sufficiently small grain sizes (10⁻² mm), the strength is sufficiently high, and the gas hydrate granule remains intact, given the above conditions. On the contrary, if the ice grain size is two orders of magnitude greater, the pressure drop destroys the granule.

To evaluate certain estimates for somewhat different values of the parameters that are closer to this study, we consider two temperature values, −10 °C and −3 °C. The gas pressure outside the granule is 0.5 MPa, while inside the granule it is equal to the equilibrium pressure for the gas hydrate, and the entire granule radius is 2.5 mm. Figure 2 shows a comparison of the maximum tensile stresses (specific expression can be found, for example, in [46,49]) in the ice crust for various external pressure values, curves (1) and (2), and ice tensile strength plotted by dashed lines (3)–(5), for three different grain sizes. These were calculated from the expression obtained on the basis of experimental data given in [48]. Horizontal axes represent the ice thickness.

From Figure 2a, it follows that the stress in the crust is higher for lower external pressure, due to increase in the pressure drop across the crust. If the ice grain size exceeds 0.1 mm (curve 5), and the external pressure is 0.5 MPa (curve 1) at temperature $-3 °\mathrm{C}$, the granule disintegrates for all crust thickness values within the modeling range. In other cases, for small thicknesses (to the left of the point of intersection of solid and dashed lines), the granule is destructed, while for larger ones, the crust withstands the load. In Figure 2b, the stress curves for two different temperatures are compared. It can be seen that the stress in the ice layer is higher for higher temperatures, while the strength is lower, which leads to an increase in the hydrate destruction range. To plot the tensile strength lines in Figure 2, we used data from [48,50].

Thus, when the external pressure is released, for any given temperature there is a certain ice thickness range, which is affected by destruction. This range becomes wider with a larger pressure drop, as well as with higher temperatures. Noticeably, granules with a sufficiently thick crust are not destroyed. The critical thickness of the ice layer can be calculated for given parameters. For example, it equals approximately 0.34 mm for the case where the external pressure drops to 0.5 MPa at −3 °C with an ice grain size of 25 µm.

When estimating the ice thickness range in which the crust is destructed under a given release of external pressure, it is important to keep in mind that at lower strain rates, the ultimate strength of the ice sample decreases, i.e., destruction occurs at lower stresses [50]. For the problem considered here, where the hydrate decomposition front moves at a rate of the order of one or a few meters per year, ice strength data for deformation rate one or two orders of magnitude lower than those considered in [48] would be more realistic. For example, the ice strength at d = 0.7 mm, −7 °C and a strain rate of 10⁻⁷ s⁻¹ is approximately half that at a rate of 10⁻⁶ s⁻¹. This means that the strength thresholds in Figure 2 will be lower than shown, and under these conditions a wider range of ice crust thicknesses will be subject to destruction.

Below, by the term hydrate decomposition, we mean the portion that dissociates under the given drop in external pressure and the selected values of other parameters.

2.2. Problem Statement

It is assumed that the frozen hydrate-bearing layer consists of ice, gas hydrate and free gas saturating small parts of pores (micropores) having supramolecular size (Figure 1). In the initial state, the micropores are not connected with each other and the gas pressure inside the pores is higher than some critical values for gas hydrates under given conditions. The pressure at the upper boundary of the reservoir is also higher than the critical one, so that the reservoir is in a metastable state due to the self-preservation effect [1,2,3].

When the pressure at the upper boundary of the formation drops below the critical value, a boundary layer with large pressure gradients emerges near the boundary due to small pore sizes [37]. Significant pressure gradients lead to the destruction of the microstructure of the boundary layer and the formation of a connected system of micropores and microchannels enabling the gas outflow. As a result, the gas pressure in the boundary layer drops to the critical pressure, which causes the hydrate depreservation and the release of gas with an abnormally high pressure, which can significantly exceed the critical value. Then, the above process of boundary destruction and hydrate decomposition extends deeper and deeper into the formation.

Thus, a downward-moving front is formed, where the microscopic pore space is being destructed, triggering hydrate decomposition into ice and gas, and the released methane is transported upward. Behind this destruction front, the permeability of the reservoir is proportional to the gas pressure gradient to some positive power [37]. Below, we analyze the qualitative and quantitative regularities of the considered processes in the isothermal approximation.

2.3. Mathematical Model Formulation

In the analysis below, we focus on the reservoir region, z > 0 (domains 1 and 2 in Figure 3). The effect of the overlying layers is taken into account by setting the pressure at the upper boundary of the reservoir. It is assumed that there is no fluid motion in the medium beneath the front, while above (behind) the front it obeys the generalized Darcy’s law [51], where the permeability is described by a power-law dependence on the pressure gradient (Figure 1). It is assumed that the process is isothermal, and the gas is thermodynamically perfect.

Thus, we start with the equations of gas motion in a low-permeable medium in the isothermal approximation following [37].

$$\frac{\partial \varphi \rho }{\partial t}+div\rho u=0$$

(1)

$$u=-\frac{k}{\mu }\nabla p$$

(2)

$$\rho =\frac{p}{R_g{T}_0},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{T}_0=const$$

(3)

Here, $\varphi$ is the reservoir porosity; $p$ is gas pressure; $\rho$ is gas density in a free state; k is reservoir permeability; $\mu$ is gas viscosity; $u$ is gas filtration velocity field; $t$ is time; $R_g$ is gas constant; $T_0$ is temperature.

In the one-dimensional case under consideration, we rewrite (1)–(3) and the boundary conditions in the form

$$\varphi \frac{\partial p}{\partial t}=\frac{\partial }{\partial z}\left(\frac{kp}{\mu }\frac{\partial p}{\partial z}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=A{\left(\frac{\partial p}{\partial z}\right)}^m,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}A=const$$

(4)

$$p(0,z)=p_0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}p(t,0)=p^0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}p(t,z_{*})=p_c,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(p_0-p_c\right){\dot{z}}_{*}=\frac{k{p}_c}{\mu }\frac{\partial p}{\partial z}|\begin{array}{l}\\ z=z_{*}\end{array}$$

(5)

Here, $p_0$ is unperturbed, or latent, gas pressure in the reservoir (pressure of the free gas after the hydrate decomposition in the absence of filtration); $p^0$ is pressure maintained at the reservoir’s upper boundary; $p_c$ is critical pressure under given conditions, triggering the hydrate decomposition when reached; the last boundary condition is the mass balance equation at the phase transition front $z_{*}(t)$.

It is essential that in (4) and (5) the permeability is proportional to the pressure gradient to a certain positive power $m$. The case $m=0$ corresponds to Darcy’s law. Let us express the parameter $A$ in terms of another constant. To achieve this, write the equality

$$\hspace{0.17em}k=A{\left(\frac{\partial p}{\partial z}\right)}^m,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}A=const,$$

(6)

or

$$\hspace{0.17em}{k}^{1/m}=A^{1/m}\frac{\partial p}{\partial z},$$

(7)

Integrating this expression over the interval $\left[0,\hspace{0.17em}\hspace{0.17em}h\right]$, where $h$ is the reservoir thickness, and raising to the power $m$, we obtain an expression for $A$:

$$\begin{array}{l}A=k_0{\left(\frac{h}{p_c-p^0}\right)}^m,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_0={\left[\frac{1}{h}{\displaystyle \underset{0}{\overset{h}{\int }}{k}^{1/m}dz}\right]}^m\\ \chi =\frac{A}{\varphi \mu }=\frac{k_0}{\varphi \mu }{\left(\frac{h}{p_c-p^0}\right)}^m\end{array}$$

(8)

Therefore,

$$\hspace{0.17em}k=k_0{\left(\frac{h}{p_c-p^0}\right)}^m{\left(\frac{\partial p}{\partial z}\right)}^m.$$

(9)

Hence, it follows that $k_0$ is the medium permeability at $m=0$ (and is expressed in permeability units). From here on, we assume that $k_0$ is a constant characterizing the filtration properties of a given porous material, and it is independent of either m or boundary conditions.

By introducing a self-similar coordinate as follows [37]:

$$p=p_{c}f(\xi ),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\xi =\frac{z}{{\left[\chi \left(t-t_0\right)\right]}^{1/\left(m+2\right)}{p}_c^{\left(m+1\right)/\left(m+2\right)}}$$

(10)

we reformulate problem (4) and (5) in a self-similar form. It yields:

$$\frac{\xi }{m+2}\frac{df}{d\xi }+\frac{d}{d\xi }\left[f{\left(\frac{df}{d\xi }\right)}^{m+1}\right]=0$$

(11)

$$f(0)=f^0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}f({\xi }_{*})=1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{{\left(\frac{df}{d\xi }\right)}^{m+1}|}_{\xi ={\xi }_{*}}=\frac{\left(f_0-1\right){\xi }_{*}}{\varphi \left(m+2\right)}$$

(12)

2.4. Problem Solution

Problems (11) and (12) are analytically unsolvable, thus we find a numerical solution. However, in the case of interest, (11) and (12) have an approximate analytical solution, which allows one to clearly see its properties.

Assuming that the first term in (11) is small, the problem can be written as

$$\frac{d}{d\xi }\left[f{\left(\frac{df}{d\xi }\right)}^{m+1}\right]=0$$

(13)

$$f(0)=f^0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}f({\xi }_{*})=1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{{\left(\frac{df}{d\xi }\right)}^{m+1}|}_{\xi ={\xi }_{*}}=\frac{\left(f_0-1\right){\xi }_{*}}{\varphi \left(m+2\right)}$$

(14)

The solution to this problem has the form

$$f(\xi )={\left[{\left(f^0\right)}^{\frac{m+2}{m+1}}+\left(1-{\left(f^0\right)}^{\frac{m+2}{m+1}}\right)\frac{\xi }{{\xi }_{*}}\right]}^{\frac{m+1}{m+2}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le \xi \le {\xi }_{*}$$

(15)

$$\frac{df}{d\xi }(\xi )=\frac{m+1}{m+2}{\left[{\left(f^0\right)}^{\frac{m+2}{m+1}}+\left(1-{\left(f^0\right)}^{\frac{m+2}{m+1}}\right)\frac{\xi }{{\xi }_{*}}\right]}^{-\frac{1}{m+2}}\left(1-{\left(f^0\right)}^{\frac{m+2}{m+1}}\right)\frac{1}{{\xi }_{*}}$$

(16)

$${\xi }_{*}={\left[\frac{\varphi \left(m+2\right)}{f_0-1}\right]}^{\frac{1}{m+2}}{\left[\frac{m+1}{m+2}\left(1-{\left(f^0\right)}^{\frac{m+2}{m+1}}\right)\right]}^{\frac{m+1}{m+2}}$$

(17)

To evaluate the accuracy of the obtained solution, we note that in (11) the expression in brackets consists of two terms of the same order. Taking (15)–(17) into account, we also obtain that the ratio of the first term in (11) to the second or third term, which we denote δ, has the order

$$\delta ~\xi {\left(\frac{df}{d\xi }\right)}^{-\left(m+1\right)}~{\xi }_{*}{}^{m+2}~{\left(f_0-1\right)}^{-1}$$

(18)

Hence, it follows that the obtained solution (15)–(17) is valid for $f_0=p_0/p_c>>1$. This is exactly the case that is of the greatest interest in the context of this study. However, as can be seen from calculations, this solution provides a reasonably good accuracy even for $f_0=1.2$, which becomes slightly worse for $f_0=1.1$ (with error being about 13% for the gas flow). Thus, the above solution can be considered as sufficiently accurate for all practically important values of $f_0$. For the case where $f_0$ is very close to unity, an approximate analytical solution was found in [52] by making use of a somewhat different method.

The gas flow through the upper boundary can be found from the expression

$$\begin{array}{l}Q=\frac{\varphi p^0}{RT}{\left(\frac{p_c{\chi }^{1/\left(m+1\right)}}{t-t_0}\right)}^{\left(m+1\right)/\left(m+2\right)}{\left(\frac{df}{d\xi }\right)}^{m+1}\\ \chi =\frac{A}{\varphi \mu }=\frac{k_0}{\varphi \mu }{\left(\frac{h}{p_c-p^0}\right)}^m\end{array}$$

(19)

This yields

$$Q=\frac{p_c}{RT}{\left(\frac{k_0{h}^m{p}_c}{\mu {\left(t-t_0\right)}^{m+1}}\right)}^{\frac{1}{m+2}}{\left[\frac{\left(m+1\right)\left(f_0-1\right)}{{\left(m+2\right)}^2}\frac{1-{\left(f^0\right)}^{\frac{m+2}{m+1}}}{{\left(1-f^0\right)}^{\frac{m}{m+1}}}\right]}^{\frac{m+1}{m+2}}$$

(20)

Thus, after combining the expressions obtained, we have

$$f(\xi )={\left[{\left(f^0\right)}^{\frac{m+2}{m+1}}+\left(1-{\left(f^0\right)}^{\frac{m+2}{m+1}}\right)\frac{\xi }{{\xi }_{*}}\right]}^{\frac{m+1}{m+2}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le \xi \le {\xi }_{*}$$

(21)

$$p(\xi )=p_{c}f(\xi ),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le \xi \le {\xi }_{*}$$

(22)

$${\xi }_{*}={\left[\frac{\varphi \left(m+2\right)}{f_0-1}\right]}^{\frac{1}{m+2}}{\left[\frac{m+1}{m+2}\left(1-{\left(f^0\right)}^{\frac{m+2}{m+1}}\right)\right]}^{\frac{m+1}{m+2}}$$

(23)

$$Q=\frac{p_c}{RT}{\left(\frac{k_0{h}^m{p}_c}{\mu {\left(t-t_0\right)}^{m+1}}\right)}^{\frac{1}{m+2}}{\left[\frac{\left(m+1\right)\left(f_0-1\right)}{{\left(m+2\right)}^2}\frac{1-{\left(f^0\right)}^{\frac{m+2}{m+1}}}{{\left(1-f^0\right)}^{\frac{m}{m+1}}}\right]}^{\frac{m+1}{m+2}}$$

(24)

$$\xi =\frac{z}{{\left[\chi \left(t-t_0\right)\right]}^{1/\left(m+2\right)}{p}_c^{\left(m+1\right)/\left(m+2\right)}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\chi =\frac{k_0}{\varphi \mu }{\left(\frac{h}{p_c-p^0}\right)}^m$$

(25)

3. Results

We solved the problem with the parameter values shown in

Table 1. Parameter values used to calculate gas flow data.

Parameter Name	Value
$\mathrm{formation} \mathrm{thickness} h$	20 m
$\mathrm{unperturbed} \mathrm{permeability} k_0$	10−23—10−18 m2
$\mathrm{gas} \mathrm{viscosity} \mu$	10−5 Pa·s
$\mathrm{temperature} T_0$	266 K
$\mathrm{gas} \mathrm{constant} \mathrm{for} \mathrm{methane} R_g$	519.65 J/kg·K
$\mathrm{critical} \mathrm{pressure} p_c$	0.85 MPa
$\mathrm{pressure} \mathrm{at} \mathrm{the} \mathrm{upper} \mathrm{boundary} \mathrm{of} \mathrm{the} \mathrm{reservoir} p^0$	0.83 MPa

.

To estimate the latent gas pressure $p_0$ in gas hydrate we assume that the composition of the specific gas hydrate is determined by the hydration number $n$ equal to the number of water molecules bound to a gas molecule. For calculations related to methane hydrate accumulations, the value $n=6$ is commonly used [3].

For the further analysis we introduce the following variable notations: ${\delta }_h$ is hydrate saturation (the proportion of hydrate in pores); ${\rho }_i$ is density of ice; ${\rho }_h$ is density of the hydrate; $s_i^0$ and $s_g^0$ are saturations with ice and gas of the pore space free from hydrates before their dissociation, $s_i^0+s_g^0=1$. Then, for $p_0$ we obtain the expression

$$p_0=\frac{\left(1-{\delta }_h\right)s_g^0{\rho }_g^0+\left(1-{\gamma }_i\right){\rho }_h{\delta }_h}{\left(1-{\delta }_h\right)s_g^0+\left(1-{\gamma }_i\frac{{\rho }_h}{{\rho }_i}\right){\delta }_h}{R}_g{T}_0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\rho }_g^0=\frac{p^0}{R_g{T}_0}$$

(26)

In this expression, ${\rho }_g^0$ and $p^0$ are the density and pressure of the gas in micropores before the dissociation; ${\gamma }_i$ is mass fraction of water (ice) molecules released during the methane hydrate dissociation.

Assuming ${\rho }_h\approx {\rho }_i$ in (26) and neglecting the density of the gas in micropores before the hydrate decomposition compared to the density after decomposition, we have

$$p_0\approx \frac{\left(1-{\gamma }_i\right){\rho }_h{\delta }_h}{\left(1-{\delta }_h\right)s_g^0+\left(1-{\gamma }_i\right){\delta }_h}{R}_g{T}_0$$

(27)

Let us specify input parameters, for example, as the following ${\delta }_h=0.8$, $s_g^0=0.1$, ${\rho }_h=800$ and ${\gamma }_i=0.87$. This means that hydrate, ice and gas occupy 80%, 18% and 2% of the pore volume, respectively. From (27), for these data we obtain the estimate $p_0\approx 90$ MPa. It follows that for the given value of the critical pressure $p_c={10}^6$ Pa, the normalized pressure $f_0=p_0/p_c$ can reach 90.

Problem (11) and (12) were also solved numerically, and the comparison of the approximate analytical and numerical solutions is given in Figure 4. Solid lines in Figure 4a show the self-similar numerical solutions for different parameter values (see the caption). The corresponding curves calculated on the basis of an approximate analytical solution are shown by dotted lines. Figure 4b shows similar curves for the solution derivative, which is proportional to the gas flow rate. It can be seen from these graphs that when the latent pressure of the gas in the hydrate $f_0=p_0/p_c$ is equal to 1.1, the analytical solution displays sufficient accuracy, and it increases further rapidly with increasing $f_0$. For smaller values, i.e., $1\le f_0<1.1$, accuracy is lower, and a numerical solution or some other kind of approximate analytical solution, for example, the one derived in [52], should be used. In this study, we are primarily interested in relatively large and large values of $f_0$, where the analytical solution provides good and very good accuracy, respectively.

One of the most important parameters of the problem is the velocity of the hydrate decomposition front ${\dot{z}}_{*}$. Figure 5a,b show the dependence of the front displacement $z_{*}$ on the latent pressure $f_0=p_0/p_c$ and the exponent $m$ (used in the expression for the permeability), respectively, a year after the start of the hydrate decomposition process. It follows from these figures that high latent pressures slow down the movement of the gas release front. The relation for this velocity reduction can be obtained from solutions (21)–(25):

$$z_{*}~{\left(f_0-1\right)}^{-1/(m+2)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{f}_0>>1$$

(28)

Increase in the exponent $m$, on the contrary, accelerates the movement of the front. However, at large times, $z_{*}$ curves for different $m$ become close to each other. Moreover, according to the analytical solution, even inversion occurs at sufficiently long times, so that larger values $z_{*}$ correspond to smaller values of the exponent $m$. However, the inversion takes place beyond the physically meaningful range of parameters, since by this time the reservoir would undergo complete decomposition, e.g., at a very low permeability coefficient $k_0={10}^{-20}$ m², and regarding other parameters relevant to Figure 5c, inversion time is about ${10}^6-{10}^7$ years.

Figure 5c shows the time dependence of the decomposition front position $z_{*}$ over a ten-year interval at $f_0=10$ for various values of the exponent $m$. According to (21–25), this dependence obeys the following expression

$$z_{*}~A(m){\left(t-t_0\right)}^{1/(m+2)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}A(m)=const$$

(29)

Although, as noted above, from this relation it follows that $z_{*}$ asymptotic for large $m$ at large time should be smaller, in reality, due to the factor $A(m)$, at larger $m$ the front moves faster up to complete decomposition of the reservoir.

Figure 5d shows the dependence of the hydrate decomposition front position on the logarithm of the relative effective permeability of the formation. Dashed lines are calculated assuming the linear Darcy’s law ($m=0$). It follows from this figure that for the considered parameters, despite the very low unperturbed permeability $k_0$, due to the dependence of the effective permeability on the gas pressure gradient, the front can travel one to two or more meters in the first year. Further, its movement slows down.

Figure 6 shows the dependence of the annual gas flow rate Q [g/m²·year] across the upper boundary of the reservoir averaged over a hundred years for different values of the parameters. From Figure 6a,b it follows that with an increase in the latent gas pressure, the flow rate increases rapidly. From the same figures it follows that with an increase in the exponent $m$, the flow rate also increases. Similar to what is observed for the phase transition (hydrate decomposition) front, at sufficiently large times, these curves converge for different $m$. Figure 6c shows the annual flow rate of gas released from the reservoir over a 100-year period; for each time moment t, the annual flow rate averaged over time interval [0, t] is plotted along the vertical axis.

Figure 6d shows the dependence of the logarithm of the annual gas flow rate averaged over 100-year period on the logarithm of the relative effective permeability for various parameters. Dashed lines indicate the corresponding behavior assuming linear Darcy’s law ($m=0$). It follows from this figure that taking into account the dependence of permeability on the pressure gradient leads to a significant increase in permeability. At the same time, for large permeability coefficients, the curves converge.

In concluding this section, we briefly discuss the connection between the proposed model and the results obtained on its basis with experimental data. The most important parameter characterizing the fluid flow in hydrocarbon production applications is the permeability to gas or oil. One of the first works assuming power-law dependence for the permeability on the pressure gradient was a gas outflow model for shale gas and oil fields [53]. This work showed a good agreement between their theoretical results and experimental data for a number of fields in the US. At the same time, the exponent $m$ for different deposits was different, varying within the interval [0.23, 1.3]. A significant monotonic dependence between effective permeability and pressure gradient during oil migration in low-permeability reservoirs based on a comparison of theory with experimental data was also confirmed in [39].

We also note that the mechanism of microstructure destruction in frozen rocks and the decomposition front formation, accompanied by intense gas release, proposed in this paper, as well as in [36], is consistent with the results of experimental studies involving agglomerate samples consisting of gas hydrate and blocks of gas-saturated ice, conducted at the Department of Geocryology, Faculty of Geology, Moscow State University [3].

From the results shown in Figure 5 it follows that for the characteristic values of the parameters, the rate of the decomposition front motion in the hydrate formation can reach several meters per year. For comparison, we refer to the field experimental data based on high-resolution seismic surveys and borehole data collected in the coastal area of the Laptev Sea in the Lena River Delta region [15], where the gas front motion rate was estimated as 5–7 m per year. Although the latter study imaged the degradation-associated gas front moving upwards, while in this paper we investigate the decomposition process propagating downwards, such a coincidence in terms of rates is apparently not accidental and is significantly determined by the permeability and changes in the physical, chemical and mechanical properties of frozen gas-containing rocks under conditions where degradation occurs.

Finally, experimental studies on the East Siberian Shelf (ESS) have shown that the average methane emission is 3 mg/(m²·day) in the main part of the ESS and reaches 13 mg/(m²·day) in the area of localized methane plumes [16]. These data are in good agreement with the gas flow rate curves shown in Figure 6.

4. Conclusions

We propose a generalized mathematical model for the decomposition of a hydrate-bearing reservoir with abnormally low permeability. The model employs the isothermal problem of the gas flow across the upper boundary of a horizontal reservoir as a result of the decomposition of a metastable gas hydrate into gas and ice. Decomposition (destruction) occurs as a result of a pressure drop below a critical value $p_c$ depending on temperature, gas hydrate texture, etc. It is assumed that the real permeability of the formation is proportional to the pressure gradient to some power $m\ge 0$. The dependences of the hydrate decomposition front velocity and the gas flow rate across the upper boundary of the reservoir depending on a number of the parameters were studied. It is shown that the latent gas pressure can significantly increase the free gas flow rate and slow down the movement of the decomposition front in a reservoir containing both gas hydrate and ice.

An approximate analytical solution of the problem in a self-similar formulation, convenient for calculations, was found. To estimate its accuracy of the solution, we also obtained a numerical solution, and showed that if latent gas pressure in the hydrate $f_0=p_0/p_c$ is equal to 1.1, the analytical solution demonstrates sufficient accuracy, and it rapidly increases with increasing $f_0$. For smaller values, i.e., when $1\le f_0<1.1$, the accuracy is lower, and one should use a numerical solution or some other approximate analytical solution, for example, the solution obtained in [42]. In this study, we are primarily interested in relatively large and large values of $f_0$, for which the analytical solution provides high and very high accuracy, respectively.