The Model of Cohesionless Sediment Blowout with an Increase in the Methane Flow Rate

Abstract

Dissociation of methane hydrates in the Arctic permafrost may lead to explosive gas emission. Methane blowout may be triggered by increasing gas flow rate at a certain depth. The mechanism of rock failure and blowout under the effect of pressurized gas was studied numerically and in laboratory experiments. The problem was formulated for the unsteady flow of compressed gas depending on the flow rate at a given depth, and pore gas pressure variations were calculated as a function of depth and time. The model parameters were chosen with reference to field data. According to the model, the input of gas to friable material at an increasing rate may lead to gas blowout and density loss propagating downward as the gas pressure exceeds the overburden pressure at some depth. The laboratory system was of the type of a Hele-Shaw cell, with small glass balls as friable material confined between two glass panels. The results of physical modeling and calculations show good agreement.

1. Introduction

Deposits of gas hydrates arouse much theoretical and practical interest. Gas hydrates are metastable ice-like solid compounds that form from water and low-molecular gas under certain thermobaric conditions [1] and consist of molecules of natural gas (mainly methane) encaged in a crystal structure of ice. As the pressure and temperature change, hydrates can dissociate, and the bound gas becomes free and is released into the atmosphere. Methane hydrates existing in shallow crust constitute a rich potential source of natural gas [2,3,4], but, on the other hand, methane is the second source of greenhouse gases, after carbon dioxide, contributing to climate change. Methane has 28 to 36 times greater contribution to the 100-year global warming potential than carbon dioxide [5,6,7,8]. Gas hydrates in subsea and onshore permafrost are prone to dissociation as the permafrost is degrading as a result of climate warming. The issues of gas hydrate accumulation and dissociation, as well as the related gas emission, have been studied quite extensively. Estimates of global methane flux are available for high-latitude tundra, lakes, and ponds [9] and for the East Siberian shelf [10]. However, the methane flux from the Arctic oceans and landmasses may be underestimated [11,12] as methane concentration in the Northern Hemisphere increases in summer due to rapid gas hydrate dissociation. Some methane flux estimates were obtained in field experiments [13,14,15,16].

The available research includes gas hydrate dissociation upon warming and decompression [17,18], 1D models of hydrate dissociation in porous media [19], and variations of porosity and permeability during hydrate formation [20,21,22]. Gas released at low rates at some depths was found to migrate toward the surface, and the gas flow may produce zones of overpressure under certain conditions. As the gas pressure reaches a critical value, the rocks above fail, and the emission may be explosive. Explosive gas emission, which produces craters in the seabed and on the land [23], was simulated for the Arctic conditions [24,25]. The model [24,25] was based on crater field observations in the Yamal Peninsula and implied decomposition of gas hydrates under changed P-T conditions. Immediately before the explosion, the liberated gas was supposed to fill a cavity confined under a permafrost lid, which was originally saturated with water.

We present numerical and laboratory experiments to study the effect of increasing gas flow rate on the process as another mechanism that may be responsible for explosive emission (blowout) of methane flowing through friable sediments.

2. Mathematical Model

The numerical experiments simulated two related successive processes: gas flow through porous sediments and gas blowout as a result of pressure increase. The models of isothermal single-phase gas flow in an isotropic porous medium are based on the mass conservation and momentum conservation (the Darcy law) principles and the pressure dependence of density, porosity, permeability, and viscosity. We used differential equations for ideal gas flow in rocks obtained in the 1940s by Leonid Leibenson [26,27,28]. The parabolic nonlinear differential equation for gas pressure is

$$\frac{k}{2\mu \left(T\right)m}{\nabla }^2{p}^2=\frac{\partial p}{\partial t}$$

(1)

The equation is valid under several assumptions: (i) the permeability $k$ and the porosity $m$ depend neither on temperature $T$ nor on pressure $p$ (i.e., the material is non-deformable); (ii) the viscosity $\mu$ is insensitive to pressure but depends on temperature; and (iii) gas flow is isothermal, and viscosity is constant at a given temperature.

Gas flowing through porous material blows out as its pressure reaches a critical value equal to overburden pressure (total weight of dry rock plus pore fluids):

$$p>p_{lit} \left(p_{lit}=\rho gh\right)$$

(2)

where $\rho$ and $h$ are, respectively, the density and thickness of overlying material (overburden), and $g$ is the gravity acceleration.

The Leibenson equation is applied to 1D single-phase gas flow in a porous medium. The modeling domain is a 1D vertical porous channel of the height $h$. The flow rate across the base of the channel ($Q$) and the pressure on its top ($p_h$) are assumed to be constant, and the modeling domain is exposed to the known constant pressure $p_0$ at $t=0$. The process is described by the system of equations:

$$\frac{\partial p}{\partial t}=\frac{\partial }{\partial y}\left(\alpha \frac{\partial p^2}{\partial y}\right), y\in \left(0,h\right), t>0, \alpha =\frac{k}{2m\mu \left(T\right)}$$

(3)

$$\begin{gathered} p\left(t=0,y\right)=p_0,\hspace{1em}\hspace{1em}y\in \left[0,h\right] \\ \frac{\partial p\left(t,y=0\right)}{\partial y}=-\frac{Q\mu \left(T\right)}{Sk},t>0 \end{gathered}$$

$$\begin{gathered} p\left(t,y=h\right)=p_h, t>0 \\ p>\rho gy \end{gathered}$$

where $S$ is the cross-section area of the gas channel.

The gas flow equation can be rewritten as

$$\frac{\partial p}{\partial t}=2\alpha \left(p\frac{{\partial }^{2}p}{\partial y^2}+{\left(\frac{\partial p}{\partial y}\right)}^2\right)$$

(4)

The space–time domain $0\le y\le h$, $0\le t\le T$ is divided by the finite-difference grid {y_i, i = 0,1,2,…,M; t_n, n = 0,1,2,…,N} with the grid spacing $\mathsf{\Delta }y$, at a time step $\mathsf{\Delta }t$. With the grid functions $P_i^n$, $P_i^{n+1}$, the system of equations becomes

$$\begin{gathered} \frac{P_i^{n+1}-P_i^n}{\mathsf{\Delta }t}=2\alpha \left(P_i^n\frac{P_{i+1}^{n+1}-2P_i^{n+1}+P_{i-1}^{n+1}}{\mathsf{\Delta }{y}^2}\right)+2\alpha {\left(\frac{P_{i+1}^n-P_{i-1}^n}{\mathsf{\Delta }y}\right)}^2, i\in \left(0,M\right), n>0 \\ {P}_i^0=p_0,\hspace{1em}\hspace{1em}i\in \left[0,M\right] \\ \frac{P_1^{n+1}-P_0^{n+1}}{\Delta y}=-\frac{Q\mu \left(T\right)}{Sk},n>0 \\ {P}_M^{n+1}=p_h, t>0 \\ {P}_i^n>\rho gy \end{gathered}$$

(5)

3. Numerical Results

The modeling started with gas flow at parameters close to the real conditions [1,11,12,13,14,15,16] (

Table 1. Model parameters.

$\mathit{h}, \mathrm{m}$	$\mathit{k}, {\mathrm{m}}^2$	$\mathrm{m}$	$\mathit{T}, \mathrm{K}$	$\mathit{Q}, {\mathrm{m}}^3/\mathit{s}$	$\mathit{S}, {\mathrm{m}}^2$	${\mathit{p}}_0, \mathrm{Mpa}$	${\mathit{p}}_{\mathit{h}}, \mathrm{Mpa}$	$\mathit{\rho }, \mathrm{kg}/{\mathrm{m}}^3$
5	$2\times {10}^{-15}$	0.2	280	${10}^{-5}$	1	0.1013	0.1013	1800

).

The depth profile of gas pressure was obtained assuming a constant pressure of 0.1013 MPa on the channel top and a constant gas flow across the base at the 5 m depth. The pressure first increases rapidly and then reaches a plateau (Figure 1).

Gas pressure and overburden pressure $P_{lit}\left(y\right)$ change with depth (Figure 2), and blowout occurs as the gas pressure exceeds the overburden pressure.

As the gas pressure reaches the critical value (black heavy line in Figure 3), blowout may occur both on the upper (Figure 3a) and lower (Figure 3b) boundaries of the modeling domain, which differ in permeability by less than one order of magnitude.

Gas pressure in the channel may never reach the overburden pressure (Figure 4a), and no blowout will occur, or it may rise to the critical value all over the channel (Figure 4b).

The time required for the gas pressure to reach the critical value depends on the gas flow rate and properties of the material at different depths (Figure 5). Correspondingly, three scenarios of rock failure (gas blowout) are possible: it may start at the base or at the top of the gas channel, or synchronously all along the channel; or it may never occur.

4. Laboratory Experiment

Laboratory modeling of gas flow was performed using a system of the type of Hele-Shaw cell (Figure 6), with two vertically oriented plane-parallel glass panels 0.7 mm apart. The space between the panels was filled with glass balls 0.2–0.3 mm in diameter, to a height of 20 cm above the cell bottom. The spacing between the panels was large enough to accommodate a double-packed layer of balls. The required density of the ball pack (1.75 g/cm³) was achieved by vibration compaction of balls poured in small batches.

The cell base (1) had a hole (2) for gas supply from a gas balloon (3), and the cell top was open to air (4). Gas pressure was measured with a transducer (5) connected to an ADC (6), gas flow rate was measured with a float-type meter (7), and changes in the material structure were recorded with a video camera (8).

The experiment included several steps:

• Gas flow through homogeneous material without change to its structure (Figure 7a);
• Formation of low-density zones propagating toward the cell bottom (Figure 7b–e). First, a few balls formed a bump, and then, as pressure was increasing further, more balls accumulated and rose above the surface, while low-density zones appeared underneath;
• Propagation of the low-density wave toward the cell bottom with loosening the whole ball mass (Figure 7f).

5. Processing of Laboratory Data

The data collected in the laboratory experiment were processed assuming the parameters of systems and materials as in

Table 2. Parameters used in data processing.

$\mathit{h}, \mathrm{m}$	$\mathit{k}, {\mathrm{m}}^2$	$\mathrm{m}$	$\mathit{\mu }, \mathrm{Pa}\times \mathrm{s}$	${\mathit{p}}_{\mathit{b}}, \mathrm{Pa}$	${\mathit{p}}_0, \mathrm{Mpa}$	${\mathit{p}}_{\mathit{h}}, \mathrm{Mpa}$	$\mathit{\rho }, \mathrm{kg}/{\mathrm{m}}^3$
0.2	$3\times {10}^{-10}$	0.4	$1.8\times {10}^{-5}$	$p_b\left(t\right)$	0.1013	0.1013	1750

.

Note that the boundary conditions included known pressure rather than gas flow rate at the model base. The pressure was constant on the top surface but increased with time at the cell bottom. Thus, the lower boundary condition was

$$p\left(t,y=0\right)=p_b\left(t\right), t>0$$

(6)

The pressure at the cell bottom was assumed to grow linearly from 0.1013 MPa until the critical level and the loss of stability. The calculations were based on empirical values of critical pressure at respective moments of time. The loss of stability began on the surface when the pressure at the cell bottom p_b(t) reached 0.1037 MPa 64 s after $t_0$ (Figure 8).

Note that the critical level of gas pressure in the laboratory experiment (point 1 in Figure 8) corresponds to the absence of failure in the numerical model. The time-dependent pressure variations in the channel (Figure 9) correlate with the overburden weight (red curve). The pressure in the channel does not exceed the overburden pressure until 64 s. At this permeability, the regime stabilizes immediately according to the new pressure at the cell bottom. That is, the curve with the new boundary condition either crosses the failure threshold immediately or the failure never occurs.

The numerical model predicts that failure starts as the critical pressure at the domain base reaches 0.1048 MPa, 95 s after $t_0$. This value differs from that of the laboratory experiment, but the failure location is the same (on the top surface).

The time required for the density loss to propagate until the cell bottom (transition from the setting of Figure 7b to that of Figure 7f) was 3 s in the numerical model and 4 s in the physical model.

In three out of four runs in the laboratory experiment, the pressure at the lower boundary changed abruptly until some constant value (0.1073 MPa), which led to failure all over the model length. This result agrees with the numerical prediction.

6. Discussion

Recently, there has been an increasing interest in studying the processes accompanying the decomposition of methanhydrates in the Arctic shelf [11,12,13,14,15]. The released methane has a significant impact on the greenhouse effect [5,6,7,8]. Observations of increasingly powerful and frequent methane breakthroughs in the zones of gas hydrate make it relevant to study the mechanism of these methane breakthroughs caused by its flow changes in the upper horizons of the Arctic permafrost. In this article, we consider the possibility of methane and sedimental weak rock blowout with an increase in the methane flow at some depth.

The suggested modeling simulates gas flow through porous material with possible blowout. In the numerical model, the gas flow is described by the equation of Leibenson. Blowout occurs as the pore gas pressure exceeds the overburden pressure. A more correct formulation of the failure (blowout) condition requires taking into account the strength of the material. However, it will not change much the resulting pressure pattern at the first approximation as the strength is included as an added constant variable which will only increase the critical pore pressure.

In the laboratory experiments, glass balls were used as a friable material, and its strength was not important. In a way, friction of balls at contact with other balls and with the cell walls may be an analog of strength. However, the calculated pressure critical for the failure was higher than that in the physical model of gas flow and blowout, and friction was neglected. The empirical and calculated critical pressure values were different possibly because the equation of Leibenson does not include changes to porosity and permeability at the pressure level above the overburden. In general, more experimental work is required, with varied properties of materials and system parameters. Such experiments are planned for the future, together with the development of the numerical model for 2D formulation and checking other criteria of gas blowout.

Note that the simple formulation of the gas flow problem and the choice of failure criterion demonstrated that the suggested mechanism of blowout triggered by gas flow rate increases as a result of gas hydrate dissociation is physically plausible and that the gas–rock system may develop instability.

It is important to discuss the differences between the presented model and the real situation. The model represents only one side of the complex process of gas emissions during the decomposition of gas hydrates. The model does not consider the process of dissociation of gas hydrates and does not consider phase transitions. The material in the experiments is a loose body, while in reality, it is not. Accounting for strength will not change the distribution of pore pressure but will allow a more accurate description of real processes. The performed study should be considered as the first approximation to the solution to the problem of gas emissions in granular media. In the next studies, it is necessary to take into account the effect of dry friction between particles. Dry friction should lead to the appearance of localization of movements in the loose mass, which are observed in reality.

7. Conclusions

We demonstrated that the suggested mechanism of methane blowout triggered by gas flow rate increase as a result of gas hydrate dissociation in the upper horizons of the Arctic permafrost is physically plausible, though more complicated models accounting for Arctic shelf conditions should be developed.

Methane flow in gas reservoirs and the conditions of gas blowout in shallow soft sediments were simulated with a problem formulation based on Leibenson’s equation.

Rock failure may follow different scenarios:

• It may start at the base of the gas channel;
• It may start at the top of the gas channel (free surface);
• It may occur synchronously all along the channel; or
• It may never occur.

The calculations show quite good agreement with the laboratory experiments, as both predict failure onset and propagation from the top surface. However, the theoretical and empirical results differ slightly in the critical pressure at the base of the modeling domain (0.1048 MPa against 0.1037 MPa) and in the time required for the failure to cover the whole domain length (3 s and 4 s). These are absolute values, and the causes of the difference can be better understood if more experimental data become available, for varied material properties and system parameters.