Subsea Methane Hydrates: Origin and Monitoring the Impacts of Global Warming

Abstract

The East Siberian Arctic shelf is the area where the largest natural gas reserves are concentrated. The formation of permafrost of the Arctic shelf during the Ice Age contributed to the emergence of a zone of stable existence of gas hydrates in the sedimentary layer, and subsequent flooding of the shelf led to its gradual degradation, the thawing of gas hydrates and the subsequent emissions of methane into the atmosphere. In the first part of the paper, we use mathematical modeling to study the processes of the formation of subsea permafrost on the Arctic shelf considering changes in the sea levels over the past 200 thousand years. Numerical simulations show the influence of climate warming up to 2200 on the degradation of subsea permafrost and the violation of the conditions for the stable existence of methane hydrates in bottom sediments using the example of the East Siberian shelf. The second part of the paper proposes a method for seismic monitoring of the state of gas hydrates based on a solution of multi-parameter inverse seismic problems. In particular, the degree of attenuation of seismic energy is one of the objective parameters for assessing the consolidation of gas hydrates: the closer they are to the beginning of decomposition, the higher the attenuation and, hence, the lower the quality factor. In this publication, we do not solve a multi-parameter inverse seismic problem for a real geological object. This would be impossible due to the lack of necessary data. Instead, we focus on substantiating the possibility of correct solutions for the problem of the reconstruction of the absorption and velocities for a viscoelastic medium in relation to the problem of monitoring the state of gas hydrate deposits. As noted in a range of publications, the thawing of gas hydrates leads to an increase in the fluid saturation of the geological medium followed by an increase in the absorption of seismic energy—that is, a decrease in the quality factor. Thus, the methods of seismic monitoring of the state of gas hydrates to predict the possibility of developing dangerous scenarios should be based on solving a multi-parameter inverse seismic problem. This publication is devoted to the presentation of this approach.

1. Introduction

The Arctic shelf is one of the largest sources of methane emissions into the atmosphere. In our opinion, one of the possible mechanisms causing this process is the decomposition of methane hydrates preserved in the permafrost zone of continental deposits of the Arctic seas during the Ice Age. A large amount of methane is concentrated in sea-bottom sediments in the form of gas hydrates [1,2], which have been formed by the accumulation of dissolved methane under certain P-T stability conditions for low temperature and high pressure in the methane hydrate stability zone (MHSZ) [2,3].

These conditions exist in marine sediments at water depths of several hundred meters and in the permafrost zone of continental sediments at high latitudes. In the ocean, gas hydrates typically appear in bottom sediments at water depths of more than 500 m; however, in the Arctic Ocean, due to low temperatures, gas hydrates can be stable at depths starting from 250 m. The shelf of the Siberian Seas (Kara, Laptev, East Siberian and Chukchi) is an area where the existence of permafrost is possibly due to the suitable P-T stability conditions. Reliable information on the feasibility of these conditions is obtained from drilling data [4,5]. The presence of permafrost in the bottom sediments provides conditions for the formation of a gas hydrate stability zone at a water depth of less than 100 m [6,7,8,9].

The gas hydrate deposits could have been formed during ice ages under subaerial conditions as a result of the lowering of the ocean level, the freezing of exposed bottom sediments and their subsequent flooding. The reserve of relict gas hydrates of the Arctic shelf is estimated at 2–1400 Gt CH4 [10,11]. These relict hydrates are sensitive to climate change, since the temperature of the bottom sediments rises after the shelf is flooded. An increase in temperature or a decrease in pressure due to changes in the sea level lead to a violation of the stability conditions for the stability of methane hydrates, which can lead to the release of large amounts of dissolved and gaseous methane [12,13]. Increased methane emissions are also expected as a result of the intensified destabilization of gas hydrates and degradation of permafrost on the Arctic shelf in the future [14].

Observational data for the years 1920–2009 in the shallow part of the shelf and the coastal zone of the Laptev Sea and the East Siberian Sea reflect a significant increase in the bottom temperature (up to 2.1 degrees), which began in the mid-1980s [15]. Several events of unprecedented warming in bottom waters thought to be on verge of freezing all year round have been observed on the central shelf in winter [16,17]. Three-dimensional numerical modelling of the Arctic Ocean state in this century has shown an increase in the temperature of the surface and bottom layers of the Siberian Arctic Seas [18,19].

One of the possible consequences of this fact is the degradation of gas hydrates associated with an increase in the temperature of the bottom water. This increase can lead to the destabilization of gas hydrate deposits, followed by the rise of a huge amount of gas. In turn, this can provoke the formation of huge gas bubbles with a diameter of several hundred meters, which can threaten both engineering structures [20] and shipping. Thus, the temperature dynamics of bottom waters near gas hydrate deposits should be monitored carefully.

The process of methane hydrates thawing is inevitably accompanied by their transition from a solid to a gaseous state, which entails an increase in the absorption of seismic energy, particularly for S-waves. Modern industrial methods of reflection seismic are successfully used to detect gas hydrates. Due to the high velocity of seismic waves in these deposits, gas hydrates appear on wave images as special objects, called Bottom Simulating Reflectors or BSRs, characterized by reversed polarity [21,22].

The BSR amplitude is determined by the acoustic impedance contrast between gas hydrates and formations in the lower part of Section [23,24]. Seismic methods are widely used in the detection of gas hydrates and have demonstrated their effectiveness in a number of cases. However, they do not allow recovery of the petrophysical parameters of methane hydrate deposits due to their extremely complex microstructure.

To overcome the limitations of seismic methods in describing the microscale structure of gas hydrate deposits, many measurements and numerical experiments have been performed to elucidate the relationships between the micro- and macroscale properties of gas hydrate accumulations. Thus, in [25], an increase in the absorption of seismic waves with an increase in the fluid saturation of rocks was established. Consequently, during the thawing of gas hydrates, which causes an increase in fluid saturation in rocks, the absorption of seismic waves will increase as well. It is this property that we propose to use to assess the state of gas hydrates.

Therefore, the study of climate change and, as a consequence, possible variations in P-T conditions, together with the localization of gas hydrate accumulations and monitoring of their properties, such as P- and S-quality factors are necessary for reliable risk prediction in relation to various parts of the shelf.

In this paper, we consider the following questions:

• Investigation of the influence of climate warming up to 2200 on the degradation of subsea permafrost and the violation of the conditions for the stable state of methane hydrates in bottom sediments on the example of the East Siberian shelf.
• Estimation of the state of permafrost and MHSZ using a numerical model of heat transfer in shelf bottom sediments and comparison of the results with standard seismic research methods based on the theory of seismic wave propagation in media with attenuation.
• Monitoring of the state of methane hydrate deposits using multi-parameter Full Waveform Inversion (FWI) of seismic data, prioritizing a reliable reconstruction of attenuation.

We use numerical modeling to assess the state of permafrost and the conditions for the stability of hydrates in bottom sediments. To determine the localization of methane hydrates, standard seismic methods are widely used; however, monitoring their physical state requires the development of fundamentally new approaches based on the solution of multi-parameter inverse seismic problems.

2. Formation of Subsea Permafrost and Methane Hydrate Stability Zones

To conduct the study, we chose an area located within the East Siberian shelf around the New Siberian Islands (see Figure 1). Our choice of this area was due to the following factors:

• The availability of data on the geological structure as well as geothermal observations up to a depth of 200 m [26].
• Increased methane emissions in this area recorded, according to both shipborne and satellite measurements [27].

To establish the geological structure and thermophysical properties of the deposits, we used data on the geological structure of the New Siberian Islands [26]; Figure 1. Their thermophysical characteristics are specified on the base of the lithological composition, geological and genetic affiliation and rock age according to the data obtained in the adjacent territories [28,29] and consistent with the data obtained in the from the borehole on the land of the New Siberian Islands [26].

2.1. Statement of the Problem and a Brief Presentation of the Method

To study the processes of formation of subsea permafrost on the Arctic shelf, a one-dimensional model of thermophysical processes in bottom sediments is used, considering phase transitions between frozen and thawed deposits [30]. The boundary condition on the surface of bottom sediments is determined by periods of transgressions-regressions, considering changes in sea level over the past 200 thousand years [26].

The sea level scenario is based on glacioeustatic curves of World Ocean level fluctuations and their transformation in accordance with regional data [31]. For New Siberia Island there are no paleotemperature data needed to construct a paleotemperature scenario. Therefore, paleotemperature isotope curves obtained from glacier cores of East Antarctica or Greenland are used, since these are the only data that reflect global climate fluctuations.

Further, they are transformed using paleotemperature reconstructions obtained in separate time intervals on the New Siberian Islands and the coastal lowlands of the north of Yakutia. The use of isotope paleotemperature curves provides for the construction of time-continuous regional curves of air and rock temperature dynamics using discrete regional data [26,30].

The intensity of the geothermal flux in all numerical simulations was considered to be equal to 60 mW/m², which corresponds to the average value of the heat flux for the specific area. To ensure that the lower boundary of the computational area does not affect the dynamics of the permafrost, we agreed the maximal depth of the geological model should be 1.5 km. To solve the Stefan problem, the method of catching the front in a node of a spatial grid was used.

Simultaneously with the simulation of the thermal state of bottom sediments, we compute the thermodynamic boundaries of the MHSZ. To determine the equilibrium pressure ($P_H$) at which methane, water or ice and hydrate can theoretically exist in phase and chemical equilibrium at a given temperature, we used the relation [3]:

$$ln\left(P_H\right)=\sum _{n=0...5}{a}_n{(T+T_D)}^n$$

where the equilibrium pressure of the hydrate $P_H$ is in MPa, the temperature T is in K, $T_D$ is the shift in the equilibrium temperature, which depends on the of salinity [30]. The coefficients $a_n$ are functions of the ground temperature and are presented in [3].

To implement the model numerically, we use the sweep method on a spatial grid with a vertical step of 0.5 m and an implicit time scheme with a step of 1 day. Modeling was performed for the shelf area of the Laptev Sea with modern depths of 20 m.

2.2. Results of the Numerical Simulation

Using the model of heat transfer in bottom sediments with changes in boundary conditions (at the upper boundary) according to the constructed paleogeographical scenario [31], changes in the thickness of the subsea permafrost and MHSZ over the past 200 thousand years (until 1950) were estimated. In numerical calculations for the period 1950–2200, the increase in the average annual temperature of bottom water was considered.

An analysis of summer observations [6] collected over the Laptev inner shelf, including the region of our interest, has shown a warming tendency in sea waters with the bottom layer temperature rising by 0.8 °C from 1984 to 2009 (0.032 °C/year). Considering further warming over Arctic shelves in the second decade of this century, we used a warming trend of 0.02 °C/year for all seasons (see Figure 2).

Preliminary modeling of the geocryological situation over a long period (200 thousand years) made it possible to establish the general patterns of its changes. The state of the permafrost depends both on the properties of the deposits that make up the geological section (the thermal and physical characteristics, moisture content and salinity) and on the combination of upper and lower boundary conditions and the thermal conductivity of the deposits.

Figure 3 shows temperature profiles in the upper 500 m of bottom sediments at different time horizons up to 2200. According to Figure 3, the maximum differences in the state of permafrost occur after the shelf flooding seven thousand years ago and the transition of continental permafrost to a subaqueous state. At the same time, the changes in the sediment temperature in the permafrost are insignificant for the last 5 thousand years.

We present the change in the permafrost thickness over the past 20 thousand years, based on the results of numerical experiments in Figure 2. During the period of the last glacial maximum (about 20 thousand years ago) and the subsequent five thousand years, the thickness of frozen deposits reached a maximum of 440 m (Figure 2). The subsequent warming and flooding of the shelf with sea water led to an increase in the temperature of bottom sediment and thawing of the subsea permafrost (Figure 2). The estimates based on the model results show that, at present, the thickness of relict permafrost in the bottom sediments of the considered shelf area is about 300 m; Figure 2 and Figure 3.

The evolution of the methane hydrate stability zone corresponds to the dynamics of the P-T conditions, i.e., the permafrost dynamics and sea level change. The maximum changes occurred during the warm period seven–six thousand years ago, which was accompanied by a sea level rise and flooding of this part of the shelf (Figure 2). The lowering of the MHSZ upper boundary roof shows approximately 60 m (Figure 2).

According to the results of numerical experiments presented in Figure 3, the mean annual surface temperature becomes positive (in °C) after 2050. By the end of the 21st century, the boundaries have an average positive temperature of more than 50 m. An increase in the temperature of bottom sediments accelerates the degradation of permafrost from the upper boundary, Figure 2. The upper boundary of the subsea permafrost deepens to a depth of 75 m by 2200.

Violation of the conditions for the formation of the MHSZ over the past six thousand years, and, consequently, the possible degradation of gas hydrates, leads to the accumulation of methane in the subpermafrost horizons of bottom sediments. The subsequent destruction of the frozen layer, which occurs in the present period and in the future, may cause an increase in the permeability of bottom sediment and the flow of methane into the atmosphere from areas of the shallow shelf.

3. Seismic Absorption

In a frozen state, gas hydrates do not practically differ from a pure elastic medium. However, as the temperature rises, they gradually change to a liquid and then to a gaseous state. Naturally, this also changes a number of their physical properties, including their ability to absorb seismic energy.

In other words, deposits of gas hydrates from an ideally elastic state pass into a viscoelastic state and then into a liquid and gaseous state. Therefore, it appears reasonable to monitor the state of gas hydrate deposits by assessing their absorption properties. Therefore, the mathematical formulation of the problem of detecting gas hydrates and monitoring their state should use the theory of seismic wave propagation in viscoelastic media [32,33]. It is the absorption variability that will provide us with information about the state of the gas hydrate deposits. Thus, we come to consider the following multi-parameter inverse problem of seismic wave propagation:

To reconstruct using multi-source–multi receiver–multi-component seismic data, we used the following elastic parameters:

• P- and S-wave propagation velocities.
• P- and S-quality factors.

Let us dwell briefly on the representation of the mathematical model of seismic wave propagation in viscoelastic media.

3.1. Mathematical Model of a Viscoelastic Medium

The attenuation of seismic energy is characterized by a quality factor:

$$Q^{-1}=\frac{1}{2\pi }\frac{\mathsf{\Delta }E}{E},$$

where E is the seismic energy per unit volume per unit time, and $\mathsf{\Delta }E$ specifies its loss per unit volume per unit time. As a result of numerous field experiments, a remarkable fact has been established [34]:

The quality factor does not depend on a frequency in the seismic frequency range for most geological media.

At present, to describe wave processes in viscoelastic media, the Generalized Linear Solid (GSLS) Model is the most widely used, which makes it possible to provide almost any value of quality factor in a given frequency range (see, e.g., [35], Chapter 5 of the book [36].). To do this, several attenuation mechanisms are used and determined by relaxation times for strains and stresses. In this paper, we use this model in the implementation of the $\tau$ method [37].

The undoubted advantage of the GSLS model lies in the possibility of introducing the quality factor independently for S- and P-waves:

• For S-shear waves

  $$\mu \left(\omega \right)={\mu }_r\left({\displaystyle 1+\frac{i{\omega }^2{\tau }^S}{L}\sum _{l=1}^L\frac{{\tau }_l}{1+i\omega {\tau }_l}}\right)={\mu }_r\left(1+{\tau }^{S}S\left(\omega \right)\right).$$

  (1)
• For P-waves

  $$\begin{array}{c}\lambda \left(\omega \right)+2\mu \left(\omega \right)=\left({\lambda }_r+2{\mu }_r\right)\left({\displaystyle 1+\frac{i{\omega }^2{\tau }^P}{L}\sum _{l=1}^L\frac{{\tau }_l}{1+i\omega {\tau }_l}}\right)\\ \\ =\left({\lambda }_r+2{\mu }_r\right)\left(1+{\tau }^{P}S\left(\omega \right)\right).\end{array}$$

  (2)

The symbol $\mathit{r}$ here means that the given value corresponds to the frequency.

Formulas (1) and (2) lead to the following representations of P- and S-quality factors when one uses L attenuation mechanisms:

$$Q_{P,S}=\frac{{\sum }_{l=1}^L\left({\displaystyle 1+{\tau }^{P,S}\frac{{\omega }^2{\tau }_l}{1+{\omega }^2{\tau }_l^2}}\right)}{{\displaystyle \omega {\tau }^{P,S}\sum _{l=1}^L\frac{{\tau }_l}{1+{\omega }^2{\tau }_l^2}}}$$

(3)

The ${\tau }_l$ parameters determine the relaxation times for stresses in the GSLS model used and have the dimensions of time. The values ${\tau }^{P,S}$ determine the values of the quality factor Q for P- and S-waves.

Thus, the construction of a mathematical model of a viscoelastic medium for a given distribution of quality factors $Q_{P,S}$ is performed in two steps:

• Calculation of the relaxation times ${\tau }_l$, ensuring the constancy of the quality factor by minimizing the functional:

  $${\tau }_l=argmin\parallel const-Q_{P,S}\left(\omega \right)\parallel$$
• Determination of the spatial distribution of the relaxation times ${\tau }^P$ and ${\tau }^S$, providing the given values of the quality factor $Q_{P,S}$ in the computational domain.

Let us now consider the wave process caused by a point source within a viscoelastic medium. In terms of temporal frequencies, we obtain the following system of equations:

$$\left\{\begin{array}{c}{\displaystyle {\omega }^2\rho u_x+\frac{\partial }{\partial x}\left({\displaystyle ((\lambda +2\mu )(1+S{\tau }^P)div\overrightarrow{u}-2\mu (1+S{\tau }^S)\frac{\partial u_z}{\partial z})}\right)+}\hfill \\ \\ {\displaystyle \frac{\partial }{\partial z}\left(\mu (1+S{\tau }^S)(\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x})\right)=F_1\left(\omega \right)G_x(\overrightarrow{x}-{\overrightarrow{x}}_0);}\hfill \\ \\ {\displaystyle \frac{\partial }{\partial z}\left(((\lambda +2\mu )(1+S{\tau }^P)div\overrightarrow{u}-2\mu (1+S{\tau }^S)\frac{\partial u_z}{\partial z})\right)+}\hfill \\ \\ {\displaystyle {\omega }^2\rho u_z+\frac{\partial }{\partial x}\left(\mu (1+S{\tau }^S)(\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x})\right)=}\hfill \\ \\ {\displaystyle F_2\left(\omega \right)G_2(\overrightarrow{x}-{\overrightarrow{x}}_0)}\hfill \end{array}\right.$$

(4)

To close the obtained system of equations, it is necessary to introduce a boundary condition on the free surface $z=0$, setting the vertical component of the stress tensor equal to zero there and also requiring the displacement vector to tend to zero at infinity in lateral directions.

3.2. Numerical Experiments on Modeling Seismic Waves in Viscoelastic Media

To evaluate how the Q-factor variability affects the seismic waves, we performed a series of numerical experiments, changing Q in the target set of layers. To achieve this, we used an elastic model that describes the behavior of the velocities and Q-factors of P- and S-waves obtained in the area marked with a red square in Figure 1 (see detailed descriptions in the works [4,26,38]).

The distribution of P-wave velocities used in the simulation is seen in Figure 4. Unfortunately, we do not know the distribution of the S-wave velocities on the same vertical profile, and thus we adopted a fairly common practice of assuming it to be equal to the speed of P-waves divided by $\sqrt{3}$.

Drilling data proved that the layer at a depth of 200 m contains gas hydrate deposits. In the initial frozen state, gas hydrates are indistinguishable in their mechanical parameters from an ideally elastic medium; however, the temperature rises, followed by the increase in fluid saturation of the medium, which leads to an increase in absorption—that is, a decrease in quality factors. It is this property of thawing gas hydrates that appears to be a prognostic sign of an approaching catastrophic methane release.

To evaluate how much Q-factor disturbances affect seismic wave fields, we performed a series of numerical experiments simulating a seismic wave field caused by a point source. In each experiment, we changed the quality factor in the layer at a depth of 200 m (see Figure 4), since it was in this layer that gas hydrates were found by drilling. For simulation, we used the previously developed finite difference schemes on the basis of the $\tau$-method, staggered grids and the software created on its basis (see [39] for more detail).

We used a Ricker pulse with a dominant frequency of 25 Hz as a probing signal radiating by sources located on the free surface. In total, we performed three experiments with different quality factors within the target layer at a depth of 200 m. Figure 5 presents the vertical components of the wave field at a distance of 500 m from the source at a point located in the middle of the target layer. Let us note that the first arrivals presented here correspond to a P-wave, which has a predominantly horizontal polarization inside the target layer.

• Seismic trace for elastic medium with no attenuation—this case corresponds to the black line. The first arrival times exactly fit those predicted by ray theory for both P- and S-waves.
• The red line corresponds to the same elastic model and the same geometry of the relative position of the source and receiver; however, in the target layer, the quality factor is 20. There are three main differences from the case of a ideal elastic medium:
  • Increase in phase speed.
  • Amplitude reduction.
  • “Spreading” of the signal, which indicates a decrease in the dominant frequency.
  • The blue curve sets the seismic trace for a Q factor of 10. The three differences noted above only become stronger.

For each of the obtained synthetic wave fields, we also estimated their global variability as the quality factor changes. The relative changes of these wave fields reached 10–13%.

Thus, based on the numerical experiments performed, we can conclude that the variability of the quality factor even in such a narrow layer introduces significant changes in the seismic wave fields. This allows us to hope for a reliable Q-factor reconstruction and, consequently, for the effective implementation of seismic monitoring of the process of thawing of gas hydrates.

4. Full Waveform Inversion for Viscoelastic Media

Let us now consider the solution of a multi-parameter inverse problem for a viscoelastic medium. Let us assume that, for a given acquisition, we know the solution of the system of Equation (4) for the frequency range (${\omega }_1,{\omega }_2$). We treat this wave field solution as the result of the action of the nonlinear operator B on the parameters of the viscoelastic medium: the Lame coefficients $\lambda ,\mu$ and the relaxation times ${\tau }^P,{\tau }^S$. Thus, the multi-parameter inverse problem for viscoelastic media is simply the inversion of the nonlinear operator B implicitly introduced by (4):

$$\overrightarrow{u}(x_r,z_r;x_s,z_s;\omega )=B[\lambda ,\mu ,{\tau }^P,{\tau }^S],{\omega }_1\le \omega \le {\omega }_2$$

(5)

Let us search for the solution of this Equation (5) by nonlinear least squares minimizing the objective functional, characterizing the discrepancy between the measured and synthetic data:

$$\mathsf{\Phi }[\lambda ,\mu ,{\tau }^P,{\tau }^S]={\parallel \overrightarrow{u}(x_r,z_r;x_s,z_s;\omega )-B[\lambda ,\mu ,{\tau }^P,{\tau }^S]\parallel }^2$$

(6)

To solve such a problem, researchers typically use local minimization techniques that reduce to computations of the gradient (5). In particular, it is worth mentioning the conjugate gradient method and its various modifications, (see [40,41]) or quasi-Newtonian methods (see Chapter 6 in the book [42]).

The kernel of gradient-type methods is the construction of an iterative process that minimizes the objective functional (6) and updates the resulting solution. Some of the simplest and most popular approaches for this are the steepest descent, conjugate gradient and quasi-Newton methods. The key idea here is the search of the next approximation along a certain direction, whose construction is based on the knowledge of the current gradient. The simplest technique here is the steepest descent in which the next value of the parameters is sought along the direction of the antigradient providing the fastest decrease in the objective functional:

$$[{\lambda }_{k+1},{\mu }_{k+1},{\tau }_{k+1}^P,{\tau }_{k+1}^S]=[{\lambda }_k,{\mu }_k,{\tau }_k^P,{\tau }_k^S]-{\alpha }_k\nabla \mathsf{\Phi }[{\lambda }_k,{\mu }_k,{\tau }_k^P,{\tau }_k^S].$$

(7)

The step length ${\alpha }_k$ is found as a result of one-dimensional minimization:

$${\alpha }_k=argmin\mathsf{\Phi }\left[\left[{\lambda }_k,{\mu }_k,{\tau }_k^P,{\tau }_k^S\right]-\alpha \nabla \mathsf{\Phi }\right].$$

(8)

However, as we have already mentioned, the steepest descent method has a number of disadvantages that significantly slow down and sometimes stop its convergence for high-dimensional problems and complex functionals. That is why various modifications are often used instead, such as various implementations of the conjugate gradient method and quasi-Newtonian methods. In this paper, we search for a minimum point using the conjugate gradient method.

Let us now describe numerical experiments on the application of full wave form inversion in a viscoelastic medium. To do this, we use the model presented in Figure 4. Let us start by obtaining the ratios for calculating the gradient.

4.1. Gradient of the Objective Functional in the Problem of Full Waveform Inversion in Viscoelastic Media

Let us start with the derivation of relations that determine the gradient of the objective functional $\mathsf{\Phi }$ in accordance with the representation (6). As a rule, the standard deviation of the discrepancy between the acquired and the synthetic data simulated for the current model is used, and thus the objective functional can be represented as a $L_2$ dot product:

$$\mathsf{\Phi }[\lambda ,\mu ,{\tau }^P,{\tau }^S]=\left(\overrightarrow{u}(x_r,z_r;x_s,z_s;\omega )-B[\lambda ,\mu ,{\tau }^P,{\tau }^S],\overrightarrow{u}(x_r,z_r,\omega )-B[\lambda ,\mu ,{\tau }^P,{\tau }^S]\right)$$

(9)

By definition, the gradient defines the linear part of the increment of the objective functional—that is:

$$\begin{array}{c}\mathsf{\Phi }[\lambda +\delta \lambda ,\mu +\delta \mu ,{\tau }^P+\delta {\tau }^P,{\tau }^S+\delta {\tau }^S]=\mathsf{\Phi }[\lambda ,\mu ,{\tau }^P,{\tau }^S]+\nabla \mathsf{\Phi }[\lambda ,\mu ,{\tau }^P,{\tau }^S]\\ \left(\delta \lambda ,\delta \mu ,\delta {\tau }^P,\delta {\tau }^S\right)+O\left(\delta {\lambda }^2,\delta {\mu }^2,{\left(\delta {\tau }^P\right)}^2,{\left(\delta {\tau }^S\right)}^2)\right)=\mathsf{\Phi }[\lambda ,\mu ,{\tau }^P,{\tau }^S]\\ {\displaystyle +2\Re \left\{\overrightarrow{u}(x_r,z_r;x_s,z_s;\omega )-B[\lambda ,\mu ,{\tau }^P,{\tau }^S],\frac{\partial B}{\partial \lambda }\delta \lambda +\frac{\partial B}{\partial \mu }\delta \mu +\frac{\partial B}{\partial {\tau }^P}+\frac{\partial B}{\partial {\tau }^S}\delta {\tau }^S\right\}}\end{array}$$

(10)

Thus, we obtain the following expression for representing the linear part of the increment of the objective functional—that is, for its gradient:

$$\begin{array}{c}{\displaystyle \nabla \mathsf{\Phi }<\delta \lambda ,\delta \mu ,\delta {\tau }^P,\delta {\tau }^S>=\left({\left(\frac{\partial B}{\partial \lambda }\right)}^{*}<\overrightarrow{u}-B>,\delta \lambda \right)+}\\ \\ {\displaystyle \left({\left(\frac{\partial B}{\partial \mu }\right)}^{*}<\overrightarrow{u}-B>,\delta \mu \right)+\left({\left(\frac{\partial B}{\partial {\tau }^P}\right)}^{*}<\overrightarrow{u}-B>,\delta {\tau }^P\right)+}\\ \\ {\displaystyle \left({\left(\frac{\partial B}{\partial {\tau }^S}\right)}^{*}<\overrightarrow{u}-B>,\delta {\tau }^S\right)}\end{array}$$

(11)

Here, $A^{*}$ denotes the adjoint to the operator A.

The obtained representation of the objective functional gradient allows us to write down the iteration process that implements the steepest descent method as follows:

$${\lambda }_{k+1}={\lambda }_k-{\alpha }_k^{\lambda }{\left(\frac{\partial B}{\partial \lambda }\right)}^{*}〈\overrightarrow{u}-B〉$$

(12)

$${\mu }_{k+1}={\mu }_k-{\alpha }_k^{\mu }{\left(\frac{\partial B}{\partial \mu }\right)}^{*}〈\overrightarrow{u}-B〉$$

(13)

$${\tau }_{k+1}^P={\tau }_k^P-{\alpha }_k^P{\left(\frac{\partial B}{\partial {\tau }^P}\right)}^{*}〈\overrightarrow{u}-B〉$$

(14)

$${\tau }_{k+1}^S={\tau }_k^S-{\alpha }_k^S{\left(\frac{\partial B}{\partial {\tau }^S}\right)}^{*}〈\overrightarrow{u}-B〉$$

(15)

4.2. Fréchet Derivative of the Operator B

We now give a representation for the formal Fréchet derivatives of the operator $B\left[\lambda ,\mu ,{\tau }^P,{\tau }^S\right]$ with respect to the parameters of the viscoelastic medium. We restrict ourselves to constructing a derivative only for the Lame parameter $\lambda$; all other derivatives are constructed according to the same scheme. In order to avoid cumbersome calculations, we write down the system of Equation (4) in general form:

$$L(\lambda ,\mu ,{\tau }^P,{\tau }^S)\overrightarrow{u}(x,z;x_s,z_s;\omega )=\overrightarrow{F}(x,z;x_s,z_s;\omega )$$

(16)

Let us now introduce the increment of the parameter $\lambda$ and the corresponding increment of the wave field $\delta \overrightarrow{u}$, and then, up to the small of the second order, we obtain the following relation connecting these increments:

$$L(\lambda ,\mu ,{\tau }^P,{\tau }^S)\delta \overrightarrow{u}(x,z;x_s,z_s;\omega )=-\frac{\partial L}{\partial \lambda }\delta \lambda \overrightarrow{u}(x,z;x_s,z_s;\omega )$$

(17)

Here, $(x_s,z_s)$ are coordinates of the point source. Let us consider now Green’s matrix $G(x,z;x_s,z_s;\omega )$, which describes, at the point $(x,z)$, a wave field from the source point of the oriented force at the point $(x_s,z_s)$. Let us rewrite Equation (17) as:

$$\left\{\begin{array}{c}{\displaystyle {\omega }^2\rho \delta u_x+\frac{\partial }{\partial x}\left({\displaystyle ((\lambda +2\mu )(1+S{\tau }^P)div\delta \overrightarrow{u}-2\mu (1+S{\tau }^S)\frac{\partial \delta u_z}{\partial z})}\right)+}\hfill \\ \\ {\displaystyle \frac{\partial }{\partial z}\left(\mu (1+S{\tau }^S)(\frac{\partial \delta u_x}{\partial z}+\frac{\partial \delta u_z}{\partial x})\right)=-\frac{\partial }{\partial x}\left(\delta \lambda (1+S{\tau }^P)div\overrightarrow{u}\right);}\hfill \\ \\ {\displaystyle \frac{\partial }{\partial z}\left(((\lambda +2\mu )(1+S{\tau }^P)div\overrightarrow{\delta u}-2\mu (1+S{\tau }^S)\frac{\partial \delta u_z}{\partial z})\right)+}\hfill \\ \\ {\displaystyle {\omega }^2\rho \delta u_z+\frac{\partial }{\partial x}\left(\mu (1+S{\tau }^S)(\frac{\partial \delta u_x}{\partial z}+\frac{\partial \delta u_z}{\partial x})\right)=}\hfill \\ \\ {\displaystyle -\frac{\partial }{\partial z}\left(\delta \lambda (1+S{\tau }^P)div\overrightarrow{u}\right).}\hfill \end{array}\right.$$

(18)

We write down the solution of the resulting system in the form:

$${\displaystyle \delta \overrightarrow{u}(x_r,z_r;x_s,z_s;\omega )=\underset{D}{\int }G(x_r,z_r;\xi ,\zeta ;\omega )\left[\begin{array}{c}\frac{\partial }{\partial x}\left[\delta \lambda (1+S{\tau }^P)div\overrightarrow{u}(\xi ,zeta;x_s,z_s;\omega )\right]\hfill \\ \\ \frac{\partial }{\partial z}\left[\delta \lambda (1+S{\tau }^P)div\overrightarrow{u}(\xi ,\zeta ;x_s,z_s;\omega )\right]\hfill \end{array}\right]d\xi d\zeta }$$

(19)

Let us suppose that the function $\delta \lambda$ is smooth and vanishes at the boundary of D:

$${\displaystyle \delta \overrightarrow{u}(x_r,z_r,x_s,z_s;\omega )=-\underset{D}{\int }\left[\begin{array}{c}{\displaystyle \frac{\partial G(x_r,z_r;\xi ,\zeta ;\omega )}{\partial x}}\hfill \\ \\ {\displaystyle \frac{\partial G(x_r,z_r;\xi ,\zeta ;\omega )}{\partial z}}\hfill \end{array}\right]\delta \lambda (1+S{\tau }^P)div\overrightarrow{u}\left(\right)d\xi d\zeta }$$

(20)

Following the definition of a derivative, we conclude, that the solution of the system of Equation (18) defines the Fréchet derivative with respect to the parameter $\lambda$. All other derivatives with respect to the parameters of the viscoelastic medium are constructed in exactly the same way. To construct adjoint operators, we find the scalar product of both parts of (20) using an arbitrary function $\varphi (x_r,z_r,x_s,z_s,\omega )$ followed by rearrangement of the integral to extract the scalar product with the function $\delta \lambda (\xi ,\eta )$.

4.3. Numerical Experiments

For numerical experiments, we use the same vertically inhomogeneous model as before (see Figure 4); however, we introduce an additional layer with P- and S-wave attenuation of $q_p$ and $Q_s$ at the depth of 200 m (see Figure 6).

When performing numerical experiments for a viscoelastic medium, we set the main task to study the coupling of parameters; therefore, we allowed for variability of both velocities and quality factors. By the connection of parameters, we mean the possibility of their false variability. Let us explain what is meant. Suppose there is a model of a viscoelastic medium with four parameters: P- and S-velocities $V_p$ and $V_s$, and P- and S-absorption $Q_p$ and $Q_s$.

Further, we assume that we know the initial approximation, which coincides with the true distribution of parameters in the medium, except for one of them, for example, $Q_p$. Then, the presence of the coupled parameters means that, in the process of solving the inverse problem, we obtain a perturbation not only for $Q_p$ but also for other parameters.

Naturally, the coupling of parameters is not allowed when monitoring the state of gas hydrates, since here it is necessary to be sure that the results obtained are true and that the obtained change in the quality factor is really associated with a change in the physical properties of the geological object under study. As can be seen from the presented results of implementation the FWI in a viscoelastic medium, our approach is free from the mentioned drawbacks.

First, let us analyze the results of Q-factor restoration for P- and S-waves as presented in the two upper panels of Figure 7. We note the confident coincidence of the true and reconstructed quality factors, both in terms of their spatial localization and absolute values as well as the absence of manifestations of their influence on the other elastic parameters. This allows us to assert that these two groups of parameters are not coupled when solving the inverse problem for viscoelastic media by the FWI method.

To conclude this section, we note the monotonic decrease in the objective functional when minimizing by the conjugate gradient method (see Figure 8). This ensures that the iterations will be not stopped before the minimum is reached.

5. Discussion

Global climate warming and anthropogenic impact associated with the intensive industrial development of the Arctic shelf led to a change in the thermobaric conditions of its bottom layers. An analysis of observations and numerical simulation results demonstrates that, in recent decades, there has been a steady increase in the temperature of bottom waters. A permanent increase in the temperature of the bottom layer leads to a change in thermobaric conditions and immersion to a greater depth of the upper boundary of the gas hydrate stability zone.

One of the undesirable consequences of this may be the thawing of permafrost, which retains gas hydrates in a solid state, and their subsequent degradation, accompanied by the release of significant volumes of methane. Therefore, the monitoring of methane hydrate deposits in order to control their condition and assess the degree of their thawing is one of the most important tasks to ensure the safe development of the shelf of the Arctic seas.

Methane hydrates located in the stability zone are close in their physical properties to ideally elastic bodies in which there is no absorption of seismic energy. However, with an increase in temperature, thawing begins, leading to a decrease in their consolidation and, consequently, to losses of seismic energy during the propagation of elastic waves. Naturally, the closer the methane hydrates are to decomposition, the higher the absorption, and the energy losses for S-waves are much higher than for P-waves.

The most universal and reliable method for solving inverse seismic problems is the Full Waveform Inversion (FWI). The results presented in the article confirm that its application makes it possible to restore, with sufficient accuracy, all the parameters of a viscoelastic medium—the velocity and quality factor of P- and S-waves. At the same time, the stability of the solution of the multi-parameter inverse seismic problem considered in this paper for a viscoelastic medium is achieved due to the availability of a priori information regarding the location of methane hydrate accumulations. To describe their localization, we performed predictive modeling focused on recovery of the location of the methane hydrate stability zone.

In this publication, we do not aim to solve a multi-parameter inverse seismic problem for a real geological object due to the lack of data necessary for this. Instead, we focus on substantiating the possibility of correctly solving the problem of reconstruction of absorption and velocities for a viscoelastic medium in relation to the problem of monitoring the state of gas hydrate deposits. As noted in the first section, the melting of gas hydrates leads to an increase in the fluid saturation of the studied volume. This, in turn, causes an increase in the absorption of seismic energy, as shown in a number of publications.

Note that one of the main difficulties in solving a multi-parameter inverse problem is the coupling of the desired parameters, which, in some cases, leads to a decrease in the accuracy of the solution. (See, for example, [43]. Acoustic multi-parameter full waveform inversion based on the wavelet method in Inverse Problems in Science and Engineering, 2021, v. 29, 2, pp. 22–247).

Thus, conducting a two-stage hybrid study allows not only localizing the stability zones of methane hydrates with acceptable accuracy but also assessing their state in order to predict possible emergency situations associated with their decomposition with methane release.

6. Conclusions

Violation of the conditions for the formation of the methane hydrate stability zone over the past six thousand years and, consequently, the possible degradation of gas hydrates at a depth of 100–165 m below the bottom has led to the accumulation of methane in frozen and cooled horizons of bottom sediments. The subsequent permafrost thaw that is currently taking place and is likely to continue in the near future may cause an increase in the permeability of bottom rocks followed by an increase in the emission of methane into the atmosphere from some areas of the shallow water shelf.

Therefore, monitoring the process of the increasing bottom temperature, which causes the thawing of permafrost, accompanied by the thawing and degradation of gas hydrates with potential massive releases of methane is an extremely important and urgent task in light of the changing climate of the Arctic caused by global warming as well as the expanding activities of geological exploration and the construction of various complex engineering structures associated with this, such as production and injection wells, drilling platforms and pipeline systems.