Effect of Methane Gas Hydrate Content of Marine Sediment on Ocean Wave-Induced Oscillatory Excess Pore Water Pressure and Geotechnical Implications

Abstract

Methane gas hydrate-bearing sediments hold substantial natural gas reserves, and to understand their potential roles in the energy sector as the next generation of energy resources, considerable research is being conducted in industry and academia. Consequently, safe and economically feasible extraction methods are being vigorously researched, as are methods designed to estimate site-specific reserves. In addition, the presence of methane gas hydrates and their dissociation have been known to impact the geotechnical properties of submarine foundation soils and slopes. In this paper, we advance research on gas hydrate-bearing sediments by theoretically studying the effect of the hydromechanical coupling process related to ocean wave hydrodynamics. In this regard, we have studied two geotechnically and theoretically relevant situations related to the oscillatory wave-induced hydromechanical coupling process. Our results show that the presence of initial methane gas pressure leads to excessively high oscillatory pore pressure, which confirms the instability of submarine slopes with methane gas hydrate accumulation originally reported in the geotechnical literature. In addition, our results show that neglecting the presence of initial methane gas pressure in gas hydrate-bearing sediments in the theoretical description of the oscillatory excess pore pressure can lead to improper geotechnical planning. Moreover, the theoretical evolution of oscillatory excess pore water pressure with depth indicates a damping trend in magnitude, leading to a stable value with depth.

1. Introduction

Thanks to several decades of research in the petroleum industry, the volumetric properties of oilfield brine are presented in the petroleum literature, with the IAPWS Formulation 1995 for the thermodynamic properties of ordinary water substances for general and scientific uses being fundamental in this regard [1]. In addition, Rowe and Chou [2] provide a correlation based on compressibility measurements designed to calculate the specific volume of brine as a function of temperature, pressure, and salinity with a reported accuracy of 0.15% over the range of temperatures from 273 to 423 °K, 0 to 25 wt% NaCl, and pressures from 0.1 to 34.3 × 10⁶ Pa. Osif [3] presented a compressibility correlation for brine using experimental data for temperatures from 200 to 270 °F and pressures from 6.9 to 138 × 10⁶ Pa, with an uncertainty in compressibility of approximately 1.5%. On the other hand, the solubility of methane in water has been reported in the petroleum literature [4,5,6], as has its solubility in brine [7,8]. Moreover, an equation of state for the methane–carbon dioxide–water system has been developed [9].

The effect of gas solubility on the volumetric properties of aqueous phases is a possibility. For instance, considering oil and water, the effect of the solubility of gas is about 10 times as large as the effect due to a change in liquid phase volume at 500 psia [10]. Moreover, the effect of natural gas exsolution on specific storage in a confined aquifer undergoing water level decline has been reported in the hydrogeology literature, where a mathematical relation has been applied to the transient simulation of groundwater discharge from a confined aquifer system [11]. Also, the effect of gas solubility on natural gas reserves and production in tight formations has been reported [12]. Therefore, in a system where the evolution of methane in the presence of brine is a possibility, an effect on the aqueous phase compressibility is imminent, and this is where gas hydrate systems become relevant.

Arctic gas hydrate deposits are estimated to hold 100–500 gigatons of carbon [13,14], which is expected to be more than 10% of the carbon in such global gas hydrate reservoirs [13], and current sophisticated climate models [15] predict a high potential of Arctic gas hydrate dissociation if bottom water temperatures increase by two degrees during the next century [14]. Thus, in the environmental literature, the contributions of upward methane seepage from gas hydrate systems to climate warming have been reported [16,17], as has its effect on excess pore water pressure in such systems [18]. Therefore, the contribution of methane evolution in seafloor sediments will likely affect the volumetric properties of pore water and any physical process related to such volumetric properties. Considering the advancement of geotechnical engineering, subsea gas infrastructure has become possible. In this regard, a submarine pipeline (also known as a marine, subsea, or offshore pipeline) is one that is laid on or below the seabed inside a trench [19,20]. Consequently, understanding the geotechnical properties of marine sediments in the presence of pore fluid with dissolved methane is critical to ensuring the stability of subsea gas/petroleum infrastructure. Already, the effect of water wave-induced excess pore water pressure on the liquefaction potential of seafloor sediments has been reported in the geotechnical literature [21,22,23], as has the effect on the stability of monopiles [24,25,26]. To effectively deal with such geotechnical problems, the theoretical basis of the evolution of excess pore water pressure under ocean waves has been extensively studied in the geotechnical literature. The oscillatory mechanism is the basic feature of wave-induced pore water pressure, which is caused by the cyclic excitation of water waves [27]. For instance, the classical diffusion model for excess pore water pressure due to wave actions proposed by Nago contains a source term that is intimately linked to the compressibility of pore water [28]. Therefore, since the solubility of methane in gas hydrate-bearing sediments impacts the volumetric properties of pore fluids, the solution or pore pressure diffusion equations containing water compressibility terms will reflect this volumetric property in addition to controlling geotechnical properties. In the research work of Nago [28] on oscillatory pore water pressure, a rigorous derivation of wave-induced excess pore water pressure in light of the classical diffusion equation was presented, and a dimensional analysis based on fundamental variables was used to study the effect of dissolved air on pore water. In the research work of Zen and Amazak [29], a classical diffusion-type equation of excess pore water pressure was also derived, assuming a one-dimensional poroelastic system, where a source due to the oscillatory wave contained the compressibility parameter. Following the derivation, a finite element approach was employed to study excess pore pressure evolution trends. In the research work of Sui et al. [30], a fully dynamic model based on Biot’s theory of poroelasticity and constitutive relations were used to investigate seabed oscillatory responses. Regarding Liao et al. [27] research work, the Reynolds-averaged Navier–Stokes equations were the theoretical foundations, and a numerical approach was used to study wave–seabed structure interactions. However, in the literature, studies dealing with wave-induced excess pore water pressure generation in marine sediments, particularly with emphasis on the effect of methane solubility in pore water and gas hydrate content on sediment compressibility and its impact on excess pore water pressure generation trends, are lacking.

In December 2018, stratigraphic tests were conducted using a gas hydrate test well Hydrate-01 drilled in the western part of the Prudhoe Bay Unit, on the Alaska continental North Slope [31], which demonstrated variable degrees of gas hydrate saturation and methane solubility in the pore fluids of marine sediments. Moreover, Taleb et al. [32] published the effect of gas hydrate on the hydromechanical properties of gas hydrate-bearing sediments. Therefore, to increase the knowledge base, in this paper, we have shown theoretically the effect of gas hydrate on the evolution of oscillatory excess pore water pressure in marine gas hydrate-bearing sediments. The theoretical findings have enabled us to discuss the geotechnical implications of such geologic systems related to submarine slope wave-induced hydromechanical process, where oscillatory excess pore pressures are a major concern. What is even more special about our research work is that it has enabled us to prove convincingly that the hydromechanical coupling phenomenon associated with oscillatory waves requires additional measures in drilling to targets associated with gas hydrate deposits in marine sediments. In this regard, the periodic increase in pore pressure due to the excess oscillatory pore water pressure requires a safety factor to be applied to drilling muds to prevent blowout.

2. Study Background

2.1. Global Distribution and Energy Potential of Natural Gas Hydrates

Gas hydrate is an ice-like crystalline form of water in combination with low-molecular-weight gases such as methane, ethane, or carbon dioxide [33], with marine sediments accounting for 99% of the global distribution, while only 1% is accounted for by the permafrost regions [34]. Figure 1 shows the global natural gas hydrates in marine and permafrost environments, with the United States leading for marine sediment deposits. Considering energy resources, gas hydrates represent greener forms of hydrocarbon resources, with a worldwide distribution and abundance that have the potential to represent yet another front in the energy industry thanks to current research in the area of extraction strategies [35]. Considering the importance of natural gas storage form, gas hydrates offer potential that makes them attractive to use for the purposes of storing and transporting energy resources as a competitor to liquefaction and condensing methods [34].

However, while gas hydrates account for a substantial energy resource, the most recent comprehensive estimates by the Global Methane Budget indicate that the annual global methane emissions are around 570 Mt, of which gas hydrates account for 40% [37]. Therefore, considering the higher global warming potential of methane relative to carbon dioxide [38], the direct emission of gas hydrates from marine sediments due to release poses imminent environmental hazards [39] (see Figure 2). In the following section, we review the importance of marine environments in light of geotechnical engineering practices related to subsea energy infrastructure developments in the energy industry.

2.2. Subsea Energy Infrastructure and Geotechnical Implications

More than a quarter of today’s hydrocarbon fuel supply is sourced from the offshore environment [40]. Consequently, marine sediments have provided a foundation for diverse geotechnical works involving the development of subsea energy infrastructure (See Figure 3). For instance, subsea manifolds have been used in the development of oil and gas fields to simplify the subsea system, which has the potential to minimize the use of subsea pipelines and risers while optimizing the fluid flow in production systems [41]. Besides subsea manifolds, marine sediments offer foundation soils for the deployment of offshore wind and tidal wave turbines (see Figure 4 and Figure 5).

Therefore, considering the role of such infrastructure as an imminent foundation load on marine sediments, the geotechnical aspects of such vital energy projects are critical to enhancing structural integrity and long-term stability [44,45,46]. In this regard, the effective stress of marine sediments depends on pore pressure, which depends partly on fluid saturation, and this is where excess pore water pressure [47] due to the gravity loads of subsea energy structures is of primary concern. However, apart from foundation loads due to structural elements, the role of wave-induced excess pore water generation has also been documented in the geotechnical literature [48,49,50,51]. Therefore, understanding the effect of such excess pore water pressure on the geotechnical properties of marine sediments is crucial to forecasting structural stability and long-term contributions to grid energy. Consequently, the roles of two types of excess pore water pressure have been identified, namely, residual and oscillatory excess pore water pressures. Hydrologically, residual excess pore water pressure is generated due to the cyclic changes in shear stress, while oscillatory components result from oscillatory motions of the surface waves [52]. In the present study, we consider the theoretical aspects of oscillatory excess pore water pressure regarding the theoretical model and the analytical solution, which is essential to predict its contribution to global excess pore water pressure generation, which also includes structural loading and residual contributions.

3. Modeling Gas Hydrate Effect on Oscillatory Excess Pore Water Pressure

3.1. Integration of Source Term into Classical Consolidation Equation

One imminent effect of pore pressure is to induce hydrodynamic flow [53], while that of coffining pressure is to induce mechanical deformation in the poroelastic structure. Therefore, in the wave-induced hydromechanical coupling process, these two separate mechanisms will contribute to the generation of excess pore water pressure that constitutes a source term. Consequently, the differential equation describing the classical one-dimensional excess pore water pressure diffusion equation considering the effect of water waves must contain the source term.

$${\displaystyle \frac{\partial \Delta u}{\partial t}}=c_v{\displaystyle \frac{{\partial }^2\Delta u}{+\partial z^2}}+\psi$$

(1)

in which $\Delta u$ is excess pore water pressure [Pa], $c_v$ coefficient of consolidation [m²s⁻¹), and $\psi$ is the source term (Pa·s⁻¹).

3.2. Application to Oscillatory Excess Pore Water Pressure Generation

Based on the mass conservation equation, Biot’s consolidation equation provides a suitable mathematical basis for the analysis of oscillatory pore pressure generation. The mass balance equation reads as follows:

$${\displaystyle \frac{\partial \in }{\partial t}}+{\displaystyle \frac{1}{Q}}{\displaystyle \frac{\partial {\Delta u}_{os}}{\partial t}}=k{\nabla }^2{\Delta u}_{os}$$

(2)

in which ${\Delta u}_{os}$ is the oscillatory excess pore water pressure [Pa], $k$ is the intrinsic permeability [m⁻²], and $Q$ is a constant [m⁻¹s].

In the hydromechanical literature, research on volumetric strain linked to excess pore water problems related to sand has confirmed that the plastic volumetric strain that has accumulated in sand, either by drained or undrained cyclic loading, dominates the increase in the liquefaction resistance of the sand. In the hydromechanical coupling process, a link exists between the rate of volumetric strain and the corresponding rate of excess pore pressure generation, which upon a minor algebraic manipulation in the present paper yields the following [54]:

$${\displaystyle \frac{{\partial \Delta u}_{os}}{\partial t}}=\xi {\nabla }^2{\Delta u}_{os}$$

(3)

where

$$\xi ={\displaystyle \frac{\left(1-{\alpha }_B{B}_{Skemp}\right)}{\left(1-2{\alpha }_{B}B\right)}}{\displaystyle \frac{kKB}{\alpha \mu }}$$

(4)

In Equation (3), ${\alpha }_B$ is Biot’s pore pressure coefficient [54], $B_{Skemp}$ is Skempton’s pore pressure coefficient, and $K$ is the undrained bulk modulus [Pa].

The oscillatory pressure on the surface of the seafloor will be transmitted to the soil mass with depth, and the only way to account for this transmission is to invoke Skempton’s pore pressure coefficient [55]. Therefore, considering the oscillatory wave-induced surface pressure on the ocean floor, a source term must be applied to Equation (3) to complete the rate of oscillatory wave-induced excess pore water pressure generation. Following Zen and Amazak [56], we write Equation (3) with a source term as follows:

$${\displaystyle \frac{\partial {\Delta u}_{os}}{\partial t}}=\xi {\nabla }^2\partial {\Delta u}_{os}+\left(\beta -1\right){\displaystyle \frac{\partial P_b}{\partial t}}$$

(5)

where

$$P_b=P_0\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{2\pi }{T_w}}t\right)$$

(6)

Here, $\beta$ is a parameter related to the compressibility of pore water, $P_0$ is the amplitude of the oscillatory wave pressure on the seafloor [Pa], and $T_w$ is the period of the wave oscillation [s].

$$\beta =1+{\displaystyle \frac{\varnothing c_w}{m_v}}$$

(7)

in which $\varnothing$ is sediment porosity, $c_w$ the compressibility of water [Pa⁻¹], and $m_v$ is the volume compressibility of water including air [56]

Hence, Equation (5) becomes the following:

$${\displaystyle \frac{\partial {\Delta u}_{os}}{\partial t}}={\xi \nabla }^2{\Delta u}_{os}+\left(\beta -1\right)P_0{\displaystyle \frac{2\pi }{T}}\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{2\pi }{T}}t\right)$$

(8)

The amplitude of the pressure on the seafloor based on the linear wave theory is given as follows [57,58]:

$$P_0={\gamma }_{sw}{A}_{amp}$$

(9)

The compressibility of water increases with dissolved gas. Therefore, assuming that the effect of the methane gas content of seafloor sediment is negligible [59], the compressibility-related parameter in Equation (7) is 1.

To account for the effect of gas hydrates within marine sediments, Equation (8) can be solved subject to two physically relevant conditions involving the presence of dissolved methane gas in contact with pore water and the presence of both dissolved and free gas. Regarding the former, the following boundary conditions hold:

$${\displaystyle \frac{\partial {\Delta u}_{os}}{\partial t}}=0atz=Hatalltimest$$

(10)

$${\Delta u}_{os}\left(0,t\right)=0$$

(11)

$${\Delta u}_{os}\left(z,0\right)=C$$

(12)

$${\displaystyle \frac{\partial {\Delta u}_{os}}{\partial z}}=0atHatalltimes$$

(13)

In Equations (10) through (13), $z$ is the depth of gas hydrate sediment [m], $H$ is the thickness of sediment [m], and $C$ is gas hydrate pressure [Pa].

In Equations (11) and (13), C is the gas pressure due to free methane gas, and H is the thickness of the gas hydrate sediment.

The solution yields the following:

$${\Delta u}_{os}\left(z,t\right)=\sum _0^{\mathrm{\infty }}{\displaystyle \frac{2\left(\left(\frac{16{\left(0.5+n\right)}^{4}C{\pi }^4}{244140625}-\frac{4{\left(0.5+n\right)}^{2}A{H}^2{\pi }^2}{15625}+B^{2}C{H}^4\right)exp\left(\frac{-{\pi }^2{\left(1+2n\right)}^2}{15625}\right)t+\left(\left(\frac{{4{\pi }^2\left(0.5+n\right)}^{4}cos(B\ast t)}{15625}\right)+B{H}^2\ast sin(B\ast t)\right)A{H}^2\right)sin\left(\frac{\pi z\left(1+2n\right)}{2H}\right)}{\left(\frac{16{\left(0.5+n\right)}^4{\pi }^4}{244140625}+B^2{H}^4\right)\left(0.5+n\right)}}$$

(14)

If the gas hydrate system does not contain free methane gas, C becomes zero, and Equation (14) becomes the following:

$${\Delta u}_{os}\left(z,t\right)=\sum _0^{\mathrm{\infty }}{\displaystyle \frac{2\left(\left(-\frac{4{\left(0.5+n\right)}^{2}A{H}^2{\pi }^2}{15625}\right)exp\left(\frac{-{\pi }^2{\left(1+2n\right)}^2}{15625}\right)t+\left(\left(\frac{{4{\pi }^2\left(0.5+n\right)}^{4}cos(B\ast t)}{15625}\right)+B{H}^2\ast sin(B\ast t)\right)A{H}^2\right)sin\left(\frac{\pi z\left(1+2n\right)}{2H}\right)}{\left(\frac{16{\left(0.5+n\right)}^4{\pi }^4}{244140625}+B^2{H}^4\right)\left(0.5+n\right)}}$$

(15)

In Equations (14) and (15), $A=\left(\beta -1\right)P_0{\displaystyle \frac{2\pi }{T}}\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{2\pi }{T}}t\right)$ and $B={\displaystyle \frac{2\pi }{T}}$.

To predict oscillatory excess pore water pressure based on the two analytical solutions, Equations (14) and (15), relevant petrophysical, hydrological, and geotechnical parameters are required. The following sections will be devoted to this task.

4. Methodology

The occurrence of gas hydrates in marine sediments originally containing pore water in a fully saturated manner and methane results in a multiphase system. Therefore, to quantify the effect of gas hydrates on the hydromechanical coupling phenomenon, porosity, hydrate saturation, and residual porosity are critical rock physics inputs. Next, the geotechnical implications of the different phases deserve to be thoroughly understood and modeled. Moreover, hydrological parameters, such as those pertaining to ocean wave dynamics, need be considered in relation to oscillatory waves. Therefore, the following sections sum up the methodology of the paper (see Figure 6).

4.1. Rock Physics Saturation Model

The residual porosity (${\varnothing }_r)$ after gas hydrate deposition can be obtained as follows [60,61]:

$${\varnothing }_r=\varnothing \left(1-S_h\right)$$

(16)

Therefore, Equation (7) becomes the following:

$$\beta =1+{\displaystyle \frac{\varnothing \left(1-S_h\right)c_w}{m_v}}$$

(17)

Within the residual pore space, we have methane and brine saturation (see Figure 7)

4.2. Geotechnical Model

The coefficient of volume compressibility ($m_v)$ is defined as follows:

$$m_v={\displaystyle \frac{\left(1+\upsilon \right)\left(1-2\upsilon \right)}{E\left(1-\upsilon \right)}}$$

(18)

The coefficient is defined as the compression of a soil layer per unit of original thickness per unit increase in effective stress on the load.

Chen et al. [63] researched the effective thermal conductivity of hydrate-bearing sediments with initial water saturation ranging from 20 to 60% and hydrate saturation ranging from 7.43 to 48.74%. In their study, the effective thermal conductivity was measured in the steady state using a high-pressure experimental setup, and the results demonstrated an increase with increasing hydrate saturation, attaining a maximum value in the saturation range of 25–36%. Therefore, we use mean values for water and hydrate saturations of 40% and 28% respectively. We calculate free gas saturation as follows:

$$S_g={1-S}_h-S_w=1-40-20=40\%$$

(19)

Using the effective medium theory, the effective Young modulus is calculated as follows [64]:

$${\displaystyle \frac{1}{K_{eff}}}={\displaystyle \frac{S_w}{K_w}}+{\displaystyle \frac{S_h}{K_h}}+{\displaystyle \frac{S_g}{K_g}}$$

(20)

Using value of bulk moduli as 7.9 GPa, 2.17 GPa, and 0.1 GPa for hydrate, formation water, and gas, respectively [62], the effective bulk modulus is calculated based on Equation (18) as 0.24 GPa.

The pore volume fraction ($\varnothing )$ is given as follows [65]:

$$\varnothing =S_g+S_w=0.40+0.2=0.60$$

(21)

The value of porosity calculated in Equation (19) falls within those reported by Wang et al. [62].

We use a medium value for Biot’s pore pressure coefficient of 0.7 [66]. The value of Skempton’s pore pressure coefficient ($B)$ was calculated using the following equation [67]:

$$B={\displaystyle \frac{\alpha }{\alpha +K\varphi \left(\frac{1}{K_{fl}}-\frac{1}{K_0}\right)}}$$

(22)

In Equation (21), $\alpha$ is Biot’s pore pressure coefficient, $K$ is the effective bulk modulus [Pa], $K_{fl}$ is the fluid modulus [Pa], and $K_0$ is the frame modulus [Pa].

Following the research work of Qadrouh et al. [68], the elastic modulus of the solid skeleton of the gas hydrate-bearing sediments is calculated using the Krief model (Krief et al.) [69], where the bulk modulus Kdry and shear modulus Gdry of the solid skeleton are expressed as follows:

$$K_{dry}=K_s{\left(1-\varnothing \right)1}^{{\displaystyle \frac{A^{*}}{\left(1-\varnothing \right)}}}$$

(23)

In Equation (23), $K_{dry}$ is the frame modulus ($K_0$) in Equation (21) [Pa], $K_s$ is the solid modulus [Pa], and $A^{*}$ is an empirical constant set to 3 [70].

Based on Equation (23) we calculate the solid modulus using effective medium theory as follows [70]:

$$K_s={\displaystyle \frac{1}{2}}\left[{\sum }_{i=1}^m{f}_i{K}_i+{\left({\sum }_{i=1}^m{\displaystyle \frac{f_i}{K_i}}\right)}^{-1}\right]$$

(24)

where m is the number of solid components, fi is the volume fraction of the i-th component, and Ki is the bulk modulus of the i-th component.

In this paper, we assume a gas hydrate-bearing marine sediment to consist solid-wise of clay and quartz, with bulk moduli given as 20.9 GPa and 36 GPa, [70] respectively, which are non-variable solid components. Assuming the volume contents of quartz and clay to be 0.8 and 0.2, respectively [70], $K_s$ is calculated as 32 GPa. Therefore, from Equation (21), and using a porosity value of 0.6 as calculated earlier, $K_{dry}=K_0$ = 16 GPa. Based on the procedure, we calculate Stefan’s pore pressure coefficient as 0.91.

In this paper, we consider gas hydrate sediments in the Bay of Fundy, Canada. Fern [71] has provided the stratigraphy of marine sediments in the Bay of Fundy (see Figure 8). In this regard, we assume glacimarine sediments with an average thickness of 30 m underlying a water depth of 150 m to host gas hydrates. For simplicity, we assume the total hydrostatic pressure at the center of the sediment to consist of 2 contributions, namely that due to the water column and that due to the weight of the sediment above the center of the sediment. We assume the submerged unit weight of the glaciomarine silt to be 11.39 kNm⁻³. The total hydrostatic pressure in this regard is 3.22 MPa.

The compresssibility of water associated with natural gas systems within the pressure range of 100 MPa is constant and close to 4.3 MPa⁻¹ [72]. Therefore, for the systems under consideration in this paper, the pressure is 4.3 × 10⁻¹⁰ Pa⁻¹.

In geophysical literature, the normalized permeability describes the ratio of the permeability of hydrate-bearing sediments to that of hydrate-free sediments. Using the bundle of parallel capillary tubes model, the permeability of the hydrate-bearing sediments is calculated as follows [73]:

$$k_n=1-S_h^2+{\displaystyle \frac{2{\left(1-S_h\right)}^2}{ln\left(S_h\right)}}$$

(25)

in which $k_n$ is the normalized permeability.

Substituting the hydrate saturation calculated earlier (0.4) into Equation (25) yields a normalized permeability of 0.107.

We assume the permeability of silty sand to be 10⁻¹⁰ ms⁻¹ (Stranne et al.) [74], which is equivalent to 1.39 × 10⁻¹⁷ m².

Thus, using Equations (24) and (25), the permeability of gas hydrate-bearing sediments can be calculated. The approach sets the intrinsic permeability to $1.48\times {10}^{-18}{\mathrm{m}}^2$.

Using the values of Biot’s pore pressure coefficient, the effective bulk modulus, and Skempton’s pore pressure coefficients determined earlier, Equation (4) can be used to calculate the epsilon parameter. Thus, assuming the viscosity of seawater at an average salinity of 35 g/kg, the viscosity of water is assumed to be 0.0014 Pa·s. The value of the parameter $\xi$ in Equation (4) is calculated using these data to be 1.47 × 10⁻⁶ m²s⁻¹.

In the calculation of the epsilon parameter, the quantity $1-2{\alpha }_{B}B$ has a negative value, which makes the value meaningless, given that it has a unit of diffusivity. Therefore, the value of the product of Biot’s pore pressure coefficient and Skempton’s pore pressure coefficient multiplied by 2 was reduced by multiplying by the porosity to arrive at a meaningful value. This approach was adopted because in the classical radial diffusion equation, the denominator of the diffusivity parameter contains a porosity multiplier [75].

The Young’s modulus E₅₀ is defined as the Young’s modulus at 50% of the maximum deviator stress. Such a Young’s modulus increases remarkably with the increase in the hydrate saturation, which is governed by the effective confining pressure. To determine the Young’s modulus, we assume an overburden pressure equal to the pressure due to the column of seawater over the sediments (see Figure 8). Assuming the density of seawater to be 1030 kgm⁻³, the pressure is 1.5 MPa. Using the Poisson [76] effect for calculating lateral stress from axial stress, the confining stress is 0.5 MPa. The calculated value enables us to determine approximately the required Young’s modulus, using Appendix A [77]. The value is deduced approximately using the plot for 0.78 MPa and assuming the hydrate concentration of 0.4 to be 120 MPa.

For most solids, the Poisson ratio is 0.25 [78]. Using Equation (16) and the value for the compressibility of water in the presence of methane gas, the value of the coefficient of volume compressibility is 6.94 × 10⁻⁹ Pa⁻¹.

4.3. Hydrological Parameters

Based on the research findings of Li et al. [79], the effect of waves generally decreases from southwest to northeast in the Bay of Fundy, and the mean significant wave height is the greatest in the outer bay and the Gulf of Maine (1–1.6 m), decreasing to 1.0–0.5 m in the mid-bay, and being reduced further to less than 0.5 m in the upper bay. Wave periods reach 6 s in the outer bay, decreasing to 4–5 s in the mid-bay and to less than 4 s in the upper bay. In another study by Amos and Asprey [80], a mean wavelength of approximately 0.3 m and a mean period of 4 s were reported. Consequently, based on the findings of the cited studies, a mean wave period of 2.50 s was used.

Using the outlined methodology, the beta parameter defined by Equation (17) is calculated to be 1.01. Using the mean wave period of 2.5 s, the $A$ and $B$ parameters in Equations (14) and (15) assume the following forms:

$$A={0.01P}_0{\displaystyle \frac{2\pi }{T}}\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{2\pi }{T}}t\right)=406\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{2\pi }{T}}t\right)$$

$$B={\displaystyle \frac{2\pi }{2.5}}=2.5$$

Equation (9) can be used to calculate the maximum pressure (P₀) due to the amplitude of the wave, using the density of seawater (1030 kgm⁻³) and a wave amplitude of 1.6 m [79]. The calculation yields a value equal to 16,150.40 Pa.

Equation (17) permits the calculation of the beta parameter with the variable compressibility of water to determine the effect of methane gas hydrate on the oscillatory excess pore water pressure defined by Equations (14) and (15).

The stiffness of uncemented sediments is determined by the effective stress regime [81], and sediments with low stiffness will experience considerable deformation, resulting in low excess pore pressure generation in the hydromechanical coupling sense. Xu [82] published experimental data on the effect of the gas hydrate volume fraction of dissociation on excess pore pressure generation as a function of effective stress (see Appendix B). In this paper, we use Appendix B to determine the effect of initial excess pore pressure C in Equation (14) on oscillatory excess pore water pressure. Accordingly, we calculate the effective stress for our system with reference to the total stress at the midpoint of the sediment of 30 m thickness (see Figure 6). In this regard, the total stress consists of the surcharge due to a water column of 150 m and the weight of sediment above the midpoint. Effective stress is then calculated as the difference between total stress and pore pressure. The approach provides effective stress as 1.533 MPa based on a submerged unit of sediment equal to 11.39 kNm⁻³. Accordingly, the initial excess pore pressure is 17 MPa based on a 0.001 volume fraction of dissociated gas hydrate. The following sections sum up the results of the application of the detailed methodology.

5. Results and Discussion

5.1. Pore Evolution

In the geotechnical literature, excess pore pressure has been predicted following the dissociation of gas hydrate-bearing sediments [83,84,85]. In the present paper, we consider the initial gas pressure due to gas hydrate dissociation to be the excess pore pressure in excess of initial hydrostatic pressure. In the poroelastic coupling process, which in the present paper is limited to oscillatory wave action, Skempton’s pore pressure coefficient is a fundamental geomechanical parameter that controls and measures the extent of pore pressure generation [86]. Moreover, the pore pressure diffusion parameter is the fundamental principal determinant in the spatiotemporal evolution of excess pore pressure. In this paper, the parameter is defined by Equation (4), which integrates both Biot’s and Skempton’s pore pressure parameters. Therefore, Equation (5) theoretically and substantially describes the oscillatory excess pore pressure, and its analytical solution must provide a better insight into excess pore pressure evolution. In the current research work, time frames for excess pore water pressure calculation were chosen to be several multiples of the mean maximum wave period of 6 s (see Section 4.3). Accordingly, Figure 9 shows the oscillatory excess pore pressure for different time frames of 3600 s, 7200s, 10,800 s, 14,400 s, and 18,000 s for the case with only dissolved methane, defined by Equation (15) for the entire thickness of the gas hydrate sediment. Figure 10 shows a similar plot over ten-meter and two-meter thicknesses, respectively, of the hydrate-bearing sediment for the sake of clarity over shallower depths for geotechnical purposes. The figures show that as time progress in hydrological processes related to wave action, the amplitudes of the oscillatory excess pore pressure decrease, reaching a maximum of 155 Pa after one hour of wave action corresponding to 3600 s. This trend is observable in Figure 9, Figure 10 and Figure 11. In all cases, oscillatory excess pore pressure evolves in a manner that indicates negative pore pressure. Moreover, excess pore pressures at all times are characterized by imminent damping of the pore pressure magnitude, approaching constant values for all times.

In the geotechnical literature, excess pore pressure due to gas hydrate dissociation has been reported [87,88,89]. Therefore, in the oscillatory wave-induced hydromechanical coupling process, substantial excess pore pressure must be possible due to excessive initial excess pore pressure caused by methane hydrate dissociation processes. Accordingly, Figure 11 shows plots of the excess pore pressure evolution for similar times as found in Figure 9 through Figure 10. Figure 12 shows the plot for the case with only dissolved methane gas a similar trend with damping amplitude while Figure 13 and Figure 14 show plots for ten-meter and two-meter sediment thicknesses, respectively, all plots being defined by Equations (14) and (15). In all cases, the amplitude of oscillatory excess pore pressure decreases with time, possessing a damping character and approaching a constant value with depth as before. However, it is noteworthy that in the case of oscillatory excess pore pressure, pore pressure far exceeds that characterized by only dissolved methane gas (see Figure 9, Figure 10 and Figure 11). In the research work of Zen and Amazak [29] on oscillatory excess pore water, their Figure 15b considers both positive and negative excess pore water pressure due to the hydromechanical coupling phenomenon, as found in Figure 9 and Figure 10 of the current research work, which supports the theoretical validity of our models. Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 shows that as the depth of the sediment increases, the oscillatory excess pore water pressure decreases in amplitude/magnitude, attaining a stable value. Figure 15 shows a plot of the excess pore pressure versus depth up to 0.75 of a meter to highlight the maximum excess pore pressure for all times. Accordingly, the figure shows decreasing maximum amplitude with time and depth of sediment.

5.2. Geotechnical Implications

Submarine landslides are some of the major geohazards in offshore settings and can affect the safety of energy infrastructure such submarine pipelines due to the imminent flow of debris [90]. In the geotechnical literature, two known processes have been cited as being geotechnically relevant in relation to stability. For instance, Wang et al. [91] discussed the stability of submarine slopes due to rapid sedimentation, while Hack et al. [92] discussed slope stability in relationship to earthquakes. In all cases, the hydromechanical coupling process is the cause of excess pore pressure buildup, which can compromise the shear strength of marine slopes that can lead to failure. A recently published article by Hussein et al. [93] shows the undrained shear strength of marine sediments versus depth for three different marine settings comprising of Blake Ridge, the Peru continental margin, and Mediterranean sites. Based on the cited reference, the undrained shear strengths for a depth of 30 are 50 kPa, 75 kPa, and 25 kPa, respectively. Stark et al. [94] also determined the undrained shear strength of a deep marine setting, reporting a shear strength value of 7 kPa for a depth of 3.5 m. Therefore, in the context of the present paper, these shear strength values are below the maximum and stable excess pore pressures predicted by Equation (14), implying imminent slope instability. Moreover, the evolution of oscillatory excess pore pressure indicates that maximum excess pore pressures for all times considered will be experienced within sediments that are one meter below the seafloor, as seen in Figure 14.

Moreover, the substantially lower values of excess pore pressure predicted by Equation (15) compared to those of Equation (14) imply that in gas hydrate-bearing sediments subjected to ocean waves, negligence of the presence of methane in the hydromechanical coupling model will lead to a gross underestimation of excess pore pressure induced by wave action, as evidenced in the plots predicted by Equation (14). Accordingly, the trends of excess pore pressure predicted by Equation (14) explain submarine slope failures reported in the geotechnical literature [95], where theoretical analyses have also shown that following the dissociation of gas hydrates, there is an increase in fluid pressure and a resulting reduction in effective stress, which, in the context of the present paper, will be aggravated by the wave-induced hydromechanical coupling process.

Negative excess pore pressures have been encountered in geotechnical engineering practices involving excavations [96]. In wave-induced hydromechanical coupling, negative excess pore pressure develops following wave troughs, which reflects excavation effects on the hydrodynamic process. In Figure 7, Figure 8 and Figure 9, the evolution of negative excess pore pressure can be accounted for in this paper by considering the absence of initial excess pore pressure due to hydrate dissociation.

5.3. Implications for Casing String Design in Drilling Engineering

To drill a hydrocarbon well to a structurally or stratigraphically defined target, casing design forms a critical aspect of the drilling program. In this regard, pore pressure and fracture pressure prognosis are essential in casing design [97]. Pore pressure or formation pressure is the pressure in the pore space due to its occupying fluid, while fracture pressure is the pressure required to fracture the formation in the direction with the least principal stress [98]. Generally, the number of casing strings required to drill to the target depth depends on the pore pressure–fracture pressure window [99], and the narrower the window, the more casing strings required and the higher the drilling cost. Therefore, the implication of the findings of the present research work is that the presence of oscillatory excess pore water pressure will periodically increase total pore pressure, and pore pressure deduced from a classical sonic log approach using a well in the vicinity of the gas hydrate deposit will not be representative unless it considers gas hydrates. Moreover, the presence of the oscillatory excess pore water pressure implies periodic increases in pore pressure with depth and a decrease in the pore pressure–fracture pressure window. Thus, mud design requires additional safety factors to cater for the increase in formation pressure due to the hydromechanical coupling phenomenon to be able to accurately predict fracture occurences [100]. Also, the obvious implication of the hydromechanical coupling phenomenon is that a periodic increase in pore pressure and decrease in the pore pressure–fracture pressure window will necessitate more casing strings, which will increase the cost of drilling.

In situ pore pressure is one of the fundamental parameters required for assessing the strength and stability of geomaterials, and its identification is essential in geotechnical planning. The conventional approach to determining excess pore water pressure requires the use of piezometers [101]. However, given the recent advances in automation/digitization in measurement while drilling (MWD) [102], the findings of the present paper highlight the potential of integrating custom-built sensors for detecting the presence of excess pore water pressure due to the hydromechanical coupling phenomenon arising from oscillatory ocean waves. Such an integral component will facilitate additional primary well control measures [103] while drilling in areas of marine sediments with gas hydrate deposits.

The source term for hydromechanical coupling (Equation (1)) depends on the compressibility of pore water. Generally, the solubility of gas in water increases its compressibility, which has the tendency to minimize the hydromechanical coupling-induced excess pore water pressure from oscillatory waves. Therefore, to minimize the oscillatory excess pore water pressure effect, basins with low geothermal gradients are required to increase the solubility of methane from gas hydrate and increase the compressibility of pore water.

In the petroleum industry, the classical approach to determining fracture pressure has to do with knowing effective stress [104,105,106]. Therefore, the presence of methane hydrate and its dissociation implies an increase in formation pressure, which reduces effective stress, leading to fracture pressure reduction. The effect will further reduce the pore pressure–fracture pressure window and increase the cost involved in the casing program. The proper estimation of fracture pressure requires the modeling of methane pressure and the hydromechanical coupling phenomenon, similar to what has been performed in the current research work. Also, the pore pressure increase due to excess pore pressure generation implies that heavier mud is required, which calls for the addition of the proper amount of hematite to water-based mud.

5.4. Summary

In this research work, we have provided theoretical derivations for the hydromechanical coupling phenomenon associated with the oscillatory wave hydromechanical coupling phenomenon in addition to discussing the geotechnical implications. We have also discussed the implications related to primary well control measures in drilling. To calculate oscillatory excess pore water pressure based on our research work, hydrological parameters pertaining to wave amplitude and frequency are required. These parameters must be readily available, spanning different time frames. Therefore, the limitations regarding the applicability of the theoretical model of the current research stems from the fact that satellite data [104] acquisition is required to gather hydrological information over an extensive period to provide reliable data for the calculation of excess pore water pressure.

6. Conclusions

The stability of submarine slopes and foundation soils is of paramount importance in subsea geotechnical engineering practice. In this regard, similar to the effect of static loading on excess pore pressure generation in foundation soils onshore, which is the direct result of the hydromechanical coupling phenomenon, the effect of wave-induced hydromechanical coupling deserve to be considered in the offshore environment. Generally, seismic wave effects are known to induce the hydromechanical coupling process, which can lead to excess pore pressure generation and a reduced foundation soil shear strength. Apart from the seismogenic effect, hydrodynamic processes related to wave action in the hydrosphere, notably in marine environments, can equally generate excess pore pressure and impact submarine slope/foundation soil instability. Moreover, the presence of gas hydrates and their dissociation are known to impact the shear strength of submarine slopes without wave action where hydrate dissociation occurs. In this paper, we theoretically studied the role of oscillatory waves in excess pore pressure generation for several time frames of wave dynamics. In this regard, we focused on two geotechnically and theoretically relevant aspects, notably the hydromechanical coupling process with initial methane gas hydrate pore pressure and the case with no initial gas pressure. Accordingly, the following points sum up the conclusion of this paper.

• The oscillatory excess pore pressures for all time frames studied in this paper exhibit a damping effect from a maximum value toward a stable value as the depth of gas hydrate-bearing sediment increases.
• The excessively high values of oscillatory excess pore pressure reported in the present paper where the presence of initial methane gas pressure is considered are an imminent precursor of the submarine slope instability and failure reported in the geotechnical literature in areas with marine wave action.
• Assuming that there is no initial methane gas hydrate, pressure leads to the underestimation of excess pore pressure that could lead to improper geotechnical engineering practices.
• The assumption that there is no initial methane gas hydrate pressure leads to negative values of oscillatory excess pore pressure during certain times of wave action, which correspond to the unloading effect associated with onshore geotechnical engineering practices related to excavation.
• The existence of oscillatory excess pore water pressure requires additional secondary well control measures that call for the maximized use of safety factors in drilling in addition to increasing the casing strings required to drill to targets depths, which increases the cost of drilling.
• The amplitude of the oscillatory excess pore water decreases with the depth of the sediment, attaining a stable value.

7. Recommendations

The findings of our research on the effect of the methane gas hydrate content of marine sediment on ocean wave-induced oscillatory excess pore water pressure and its geotechnical implications underscore the need to consider not only the impact on geomechanical properties but also the impact on drilling and casing design. Therefore, to increase efficiency, investments in equipment related to satellite data acquisition in marine environments must be considered a priority. The approach will provide quality satellite data on hydrological parameters pertaining to wave height and frequency critical for assessing the hydromechanical coupling effects.