An All-At-Once Newton Strategy for Marine Methane Hydrate Reservoir Models

Abstract

The migration of methane through the gas hydrate stability zone (GHSZ) in the marine subsurface is characterized by highly dynamic reactive transport processes coupled to thermodynamic phase transitions between solid gas hydrates, free methane gas, and dissolved methane in the aqueous phase. The marine subsurface is essentially a water-saturated porous medium where the thermodynamic instability of the hydrate phase can cause free gas pockets to appear and disappear locally, causing the model to degenerate. This poses serious convergence issues for the general-purpose nonlinear solvers (e.g., standard Newton), and often leads to extremely small time-step sizes. The convergence problem is particularly severe when the rate of hydrate phase change is much lower than the rate of gas dissolution. In order to overcome this numerical challenge, we have developed an all-at-once Newton scheme tailored to our gas hydrate model, which can handle rate-based hydrate phase change coupled with equilibrium gas dissolution in a mathematically consistent and robust manner.

1. Introduction

Methane hydrates constitute a dominant organic carbon pool in the Earth system, and an important intermediate “capacitor” in the global methane budget. Gas hydrates (gas hydrates and methane hydrates are used interchangeably, as in the settings of our interest gas hydrates are predominantly composed of methane; the contribution of other gases is negligible) are predominantly formed from biogenic methane generated by the microbial degradation of organic matter (methanogenesis) in the deep biosphere. This methane migrates upwards as free gas or methane-rich porewater by advection. This fluid flow is caused by non-steady-state sediment compaction (passive margins), the compaction of oceanic sediments during subduction (active margins), and the dewatering of minerals at elevated temperatures (passive + active margins). Over geological time, the hydrates accumulate close to the bottom simulation reflector (BSR, the lower stability limit of gas hydrates) because the methane flux from below leads to hydrate formation in the gas hydrate stability zone (GHSZ). However, the ongoing sedimentation tends to bury the hydrates below the GHSZ where the hydrates dissociate, and the released methane gas migrates back into the GHSZ to re-form the hydrates. Towards the seafloor, the hydrates dissolve due to undersaturation of porewaters as a consequence of the anaerobic oxidation of methane (AOM). Some methane gas by-passes the GHSZ and AOM zone if the upward flow is larger than the reaction rates. This methane fuels rich cold seep ecosystems [1]. One of our main motivations for modeling methane hydrate geosystems is to understand the role of gas migration through the GHSZ in the natural carbon cycle [2]. Methane hydrates are also seen as an attractive future energy resource. It is estimated that the total carbon content of gas hydrates is possibly larger than the combined carbon content of all other fossil fuels [3,4,5]. However, there are a number of serious environmental risks associated with the exploitation of gas hydrate reservoirs for the purpose of gas production, such as methane leakage and related global warming, ecological hazards, and geomechanical and structural instabilities [6,7]. Therefore, we also apply our gas hydrate models for feasibility analysis and risk assessment of various gas production scenarios. It is therefore evident that we require our gas hydrate models to be robust and stable across a range of spatial and temporal scales.

One of the main challenges in modeling methane hydrate geosystems comes from the complex phase transitions which cause phases to appear and disappear locally. For example, when methane hydrates dissociate due to changes in the local thermodynamic state (i.e., temperature, pressure, and salinity), they release methane and water. Methane is released as microscopic gas bubbles which, depending on the local solubility limit for methane dissolution, either collapse into the water phase or coalesce, leading to the appearance a free gas phase. Conversely, when methane hydrates form, given the right temperature, pressure, and salinity conditions, they consume methane, which may lead to the disappearance of the free gas phase. In multiphase multicomponent models, the appearance or disappearance of phases causes the model to degenerate, such that the Jacobian becomes locally indeterminate [8,9]. For methane hydrate models in particular, additional numerical challenges arise. Firstly, the hydrate and gaseous-aqueous phase transitions are strongly coupled through nonlinear mass and thermal source and sink terms which are highly sensitive to the local thermodynamic state. Even a small change in temperature, pressure, or salinity (e.g., due to climate change, sedimentation, etc.) can locally destabilize hydrates and result in a free gas pocket with lower salinity (due to the release of fresh water bound to the hydrates), lower temperature (due to exothermic hydrate dissociation), and elevated pore pressure (due to the release of a large volume of gas). These temperature, pressure, and salinity changes may in turn push the hydrates back into the stability zone. It is therefore very difficult to determine the local phase states. Secondly, for the problems at geological scales, the rate of hydrate phase change is often several orders of magnitude lower compared to the rate of methane dissolution. This difference in the rates of the phase changes makes it even harder to determine the local phase states consistently.

A variety of numerical methods have been developed to handle the phase transitions in porous media models, such as primary variable switching (PVS) schemes [10,11], the method of negative saturations [12], the method of persistent variables [13,14], and semi-smooth Newton approaches [15,16,17,18]. In the most widely used gas hydrate reservoir simulators (e.g., TOUGH-Hydrate [19], HYRES-C [20,21], STOMP-HYD [22], HRS [23,24], including our own models [25]), the gaseous–aqueous phase transitions are handled using the PVS schemes, where the choice of the primary variables is adapted locally to the phase state. This method is conceptually relatively simple and works well for oil and gas reservoir models. However, in gas hydrate models, the phase states tend to switch back and forth rapidly due to the strong coupling and nonlinearities, which often leads to spurious oscillations and a drastic reduction in time-step size, and in the extreme case, even to a breakdown of the numerical algorithm. Some gas hydrate models (e.g., [26,27], based on Code_Bright) use the method of persistent variables. This method is more stable than PVS in general, but for gas hydrate models it can be somewhat limiting because the capillary pressure relationships for hydrate-bearing sediments are very complex functions of not only the water saturation, but also the hydrate saturation (and possibly, the porosity as well as ice saturation), and therefore some parameterizations for capillary pressure may not be easy to invert, or in extreme cases, capillary pressure may not be invertible at all (e.g., experimental measure fields).

In this article, we present a robust implicit semi-smooth Newton scheme based on an NCP approach for handling the phase transitions in our methane hydrate model which includes rate-based hydrate phase change coupled with the vapor–liquid equilibrium across the gas–water interface. The advantage of an NCP approach is that it ensures that the primary variables of our mathematical model remain the same throughout the simulation, and that the constraints are realized in a variationally consistent manner, resulting in a more robust numerical scheme. As a general outline, we cast the inequality constraints arising from the vapor–liquid-equilibrium (VLE) assumption (e.g., [28]) for the ${CH}_4-H_{2}O$ system into a set of complementarity conditions which lead to the mathematical structure of a variational inequality (e.g., [29,30]). We reformulate the complementarity conditions as a set of non-differentiable but semi-smooth functions which are solved together with the governing PDEs of the methane hydrate model fully implicitly using a semi-smooth Newton method (e.g., [31] and the references therein). We implement our semi-smooth Newton method using an active-set strategy (e.g., [32,33]), where the Jacobian is uniquely determined based on the local phase states which are partitioned into active/inactive sets using the semi-smooth NCP functions.

2. Mathematical Model

Here, we introduce our complete mathematical model for methane hydrates in the marine subsurface (extended from [25]). The model is founded on the theory of porous media on the continuum scale. The representative elementary volume (REV) underlying the model is shown in Figure 1.

The model considers two fluid phases: gaseous and aqueous; and two solid phases: porous granular material (sediment) and methane hydrate. The sediment phase is referred to as the primary matrix, the hydrate + sediment as the composite matrix, and the fluid and the hydrate phases together as the pore-filling phases The phases are identified with the subscripts g, w, s, and h, respectively. The sediment is assumed to be perfectly rigid (geomechanics is ignored within the scope of this work). The fluid phases are mobile, while the hydrate phase is assumed immobile. Methane hydrate phase change is modeled as a kinetic reaction which is strongly dependent on the local thermodynamic state of the system. The hydrate phase is assumed to contain only pure methane hydrate. Gas adsorption/desorption on the surface of hydrates is not considered. In marine settings, water salinity has a strong influence on the thermodynamics and phase transitions. Therefore, the model also considers the transport of dissolved salts. The miscibility of methane is also accounted for. Therefore, in our model the gaseous phase is comprised of two components: methane and water. Meanwhile, the aqueous phase is comprised of three components: methane, water, and salts. The components are identified by the superscripts ${CH}_4$, $H_{2}O$, and c, for methane, water, and salts, respectively. The model also accounts for the thermal effects, especially the volumetric heat generation due to hydrate phase change, but assumes a local thermal equilibrium within an REV, such that a single average temperature can be defined over an REV.

In the subsequent text, the following notation is used: Subscript “$\alpha =g,w$” denotes the fluid phases, “$\beta =g,w,h$” denotes the pore-filling phases, and the superscript “$\kappa =C{H}_4,H_{2}O,c$” denotes the components. Phase saturations are denoted with $S_{\beta }$, and mole fractions of each component $\kappa$ in each fluid phase $\alpha$ are denoted with ${\chi }_{\alpha }^{\kappa }$. Note that ${\chi }_g^c=0$. Fluid phase pressures are denoted with $P_{\alpha }$, temperature is denoted with T, and porosity is denoted with $\varphi$ ($\varphi$ refers to the total porosity, which indicates the void spaces within the primary sediment matrix. This is different from apparent porosity ${\varphi }_{eff}$, which indicates the actual void spaces available for the fluid flow. See Figure 1 for more details).

2.1. Governing Equations

2.1.1. Mass, Momentum, and Energy Conservation

The transport processes characterizing the gas production from sub-surface methane hydrate reservoir are described by invoking the conservation laws for mass, momentum, and energy described for the macroscale properties of the porous medium derived using local volume averaging principles [34,35,36]. Mass balance is considered component-wise for each $\kappa =C{H}_4,H_{2}O$,

$$\begin{array}{cc}\hfill & \sum _{\alpha }{\partial }_t\left(\varphi {\rho }_{\alpha }{\chi }_{\alpha }^{\kappa }{S}_{\alpha }\right)+\sum _{\alpha }\nabla ·\left({\rho }_{\alpha }{\chi }_{\alpha }^{\kappa }{\mathbf{v}}_{\alpha }\right)=\sum _{\alpha }\nabla ·\left(\varphi S_{\alpha }{\mathbf{J}}_{\alpha }^{\kappa }\right)+{\dot{g}}^{\kappa },\hfill \end{array}$$

(1)

where ${\mathbf{v}}_{\alpha }$ denotes the velocity of the fluid phase $\alpha$ relative to the primary sediment matrix, and ${\mathbf{J}}_{\alpha }^{\kappa }$ denotes the diffusive flux of the component $\kappa$ in the phase $\alpha$. Mass balance for the hydrate phase is given by

$$\begin{array}{c}\hfill {\partial }_t\left(\varphi {\rho }_h{S}_h\right)={\dot{g}}_h.\end{array}$$

(2)

The terms ${\dot{g}}^{C{H}_4}$, ${\dot{g}}^{H_{2}O}$, and ${\dot{g}}_h$ denote the volumetric source terms resulting from the hydrate phase change, s.t., ${\dot{g}}^{C{H}_4}+{\dot{g}}^{H_{2}O}+{\dot{g}}^h=0$. Mass balance for the dissolved salt is given by

$$\begin{array}{c}\hfill {\partial }_t\left(\varphi {\rho }_w{S}_w{\chi }_w^c\right)+\nabla ·\left({\rho }_w{\chi }_w^c{\mathbf{v}}_w\right)=\nabla ·\left(\varphi S_w{\mathbf{J}}_w^c\right),\end{array}$$

(3)

where ${\mathbf{J}}_w^c$ is the Fickian diffusion flux of salt in the aqueous phase. Momentum balance for the fluid phases can be reduced to Darcy’s Law (e.g., [28]),

$$\begin{array}{c}\hfill {\mathbf{v}}_{\alpha }=-K\frac{k_{r\alpha }}{{\mu }_{\alpha }}\left(\nabla P_{\alpha }-{\rho }_{\alpha }\mathbf{g}\right),\end{array}$$

(4)

where K denotes the intrinsic permeability of the composite matrix, $k_{r\alpha }$ denotes the relative permeabilities, and ${\mu }_{\alpha }$ the dynamic viscosities. The primary matrix is assumed rigid and the hydrate phase is assumed immobile. Therefore, the momentum of the solid phases is always conserved. For describing the energy conservation in the porous medium, one energy balance equation is sufficient since local thermal equilibrium has been assumed (e.g., [28]). The energy balance is given by

$$\begin{array}{cc}\hfill & {\partial }_t\left[\left(1-\varphi \right){\rho }_s{u}_s+\sum _{\beta }\left(\varphi S_{\beta }{\rho }_{\beta }{u}_{\beta }\right)\right]+\sum _{\alpha }\nabla ·\left({\rho }_{\alpha }{\mathbf{v}}_{\alpha }{h}_{\alpha }\right)=\nabla ·k_{eff}^{th}\nabla T+{\dot{Q}}_h,\hfill \end{array}$$

(5)

where ${\dot{Q}}_h$ denotes the heat of hydrate phase change. The term $h_{\alpha }$ is the specific enthalpy of fluid phase $\alpha$, $u_{\gamma }$ is the specific internal energy of any phase $\gamma =g,w,h,s$, and $k_{eff}^{th}$ is the effective (or lumped) thermal conductivity, s.t.

$$\begin{array}{c}\hfill h_{\alpha }={\int }_{T_{ref}}^{T}C{p}_{\alpha }dT,u_{\gamma }={\int }_{T_{ref}}^{T}C{v}_{\gamma }dT,\mathrm{and},k_{eff}^{th}=\left(1-\varphi \right)k_s^{th}+\sum _{\beta }\left(\varphi S_{\beta }{k}_{\beta }^{th}\right).\end{array}$$

2.1.2. Closure Relationships

The phase saturations and the phase pressures are not independent. The saturations of the pore-filling phases are related through the summation condition,

$$\begin{array}{c}\hfill \sum _{\beta }{S}_{\beta }=1.\end{array}$$

(6)

The pressures of the fluid phases are related through a capillary pressure $P_c$ as

$$\begin{array}{c}\hfill P_g-P_w=P_c(S_w,S_h,\varphi ).\end{array}$$

(7)

This pressure difference occurs across the gaseous and aqueous phase interface due to balancing of cohesive forces within the liquid and the adhesive forces between the liquid and soil matrix. The parametrization used for approximating $P_c$ is discussed in Section 2.2.4.

2.2. Constitutive Relations

The nine governing Equations (1)–(7) consist of the following 25 unknowns:

$$S_{\beta },{\chi }_{\alpha }^{\kappa },P_{\alpha },P_c,T,{\mathbf{v}}_{\alpha },{\mathbf{J}}_{\alpha }^{\kappa },{\dot{g}}^{C{H}_4},{\dot{g}}^{H_{2}O},{\dot{g}}_h,{\dot{Q}}_h.$$

To close the model, we define 16 additional constitutive relationships in this section for the unknowns ${\chi }_{\alpha }^{\kappa }$, ${\mathbf{J}}_{\alpha }^{\kappa }$, $P_c$, ${\dot{g}}^{C{H}_4}$, ${\dot{g}}^{H_{2}O}$, ${\dot{g}}_h$, and ${\dot{Q}}_h$. Some other properties which are important for modelling hydrate geosystems are also discussed.

2.2.1. Vapor–Liquid Equilibrium

Methane and water components are assumed to exist in a state of vapor–liquid equilibrium (VLE), and Henry’s law and the Raoult’s law are assumed to be valid:

$$\begin{array}{cc}\hfill \mathrm{Henry}’\mathrm{s}\mathrm{law},z^{C{H}_4}{\chi }_g^{C{H}_4}{P}_g& =H_w^{C{H}_4}{\chi }_w^{C{H}_4},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \mathrm{Raoult}’\mathrm{s}\mathrm{law},{\chi }_g^{H_{2}O}{P}_g& =P_{sat}^{H_{2}O}{\chi }_w^{H_{2}O},\hfill \end{array}$$

(9)

where $z^{C{H}_4}$ is the methane gas compressibility, $H_w^{C{H}_4}$ is the pressure-corrected Henry’s law solubility constant for methane dissolution in water, and $P_{sat}^{H_{2}O}$ is the saturation vapor pressure for water in contact with methane gas. In addition to relationships (8) and (9), we observe that within each phase $\alpha$, the sum of the constituent mole fractions is bounded from above by one, and the equality holds only if the phase is present:

$$\begin{array}{c}\hfill \sum _{\kappa }{\chi }_{\alpha }^{\kappa }\le 1\forall \alpha \mathrm{and}\sum _{\kappa }{\chi }_{\alpha }^{\kappa }=1\mathrm{iff}{S}_{\alpha }>0.\end{array}$$

(10)

We can cast the conditions in (10) as a set of Kharush–Kuhn–Tucker complementarity conditions [37] as

$$\begin{array}{c}\hfill 1-\sum _{\kappa }{\chi }_{\alpha }^{\kappa }\ge 1,S_{\alpha }\ge 0,S_{\alpha }\left(1-\sum _{\kappa }{\chi }_{\alpha }^{\kappa }\right)=0,\forall \alpha .\end{array}$$

(11)

2.2.2. Diffusive Mass Flux

The diffusive solute flux through the composite sediment matrix is evaluated using Fick’s Law (e.g., [28]),

$$\begin{array}{cc}\hfill & {\mathbf{J}}_{\alpha }^{\kappa }=-\tau D_{\alpha }^{\kappa }\left({\rho }_{\alpha }\nabla {\chi }_{\alpha }^{\kappa }\right),\hfill \end{array}$$

(12)

where $\tau$ denotes the tortuosity of the composite sediment matrix, and $D_{\alpha }^{\kappa }$ are the molecular diffusion coefficients for components $\kappa$ through fluid phases $\alpha$. Additionally, the summation conditions ${\displaystyle \sum _{\kappa }}{\mathbf{J}}_{\alpha }^{\kappa }=0$ hold for all phases $\alpha$. Note that ${\mathbf{J}}_g^c=0$ since ${\chi }_g^c=0$.

2.2.3. Hydrate Phase Change Kinetics

When solid methane hydrates are warmed or depressurized, they decompose into methane gas and liquid water, and vice versa. This chemical reaction is expressed as $C{H}_4·N_h{H}_{2}O\rightleftharpoons C{H}_4+N_h{H}_{2}O$, where $N_h$ gives the stoichiometry of water molecules per molecule of gas (i.e., the hydration number). The rate of this reaction is modeled by the Kim–Bishnoi kinetic model [38], where the rate of methane and water generated as a result of hydrate phase change are evaluated as

$$\begin{array}{cc}\hfill & {\dot{g}}^{C{H}_4}=k^r{M}^{CH4}{A}_{rs}\left(P_e-P_g\right),\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & {\dot{g}}^{H_{2}O}={\dot{g}}^{C{H}_4}{N}_h\frac{M^{H_{2}O}}{M^{C{H}_4}},\hfill \end{array}$$

(14)

where $P_e$ is the equilibrium pressure for the methane hydrate, $k^r$ is the kinetic rate constant, and $A_{rs}$ is the specific reaction surface area. Note that if $\left(P_e-P_g\right)>0$, the hydrate becomes unstable. In this case $k^r$ refers to the dissociation rate constant. Alternatively, if $\left(P_e-P_g\right)<0$, the hydrate becomes stable. In this case, $k^r$ refers to the formation rate constant. $M^{\kappa }$ denotes the molar weights, and for methane hydrate, $M_h=M^{C{H}_4}+N_h{M}^{H_{2}O}$. Additionally, the following condition holds:

$$\begin{array}{c}\hfill {\dot{g}}^{C{H}_4}+{\dot{g}}^{H_{2}O}+{\dot{g}}_h=0.\end{array}$$

(15)

For the hydrate phase change, the following constraints are considered: Hydrate dissociation can occur only when hydrate is available and the gas phase pressure is lower than the hydrate equilibrium pressure, and conversely, hydrate formation can occur only when both water and gaseous methane are available and the gas pressure is higher than the hydrate equilibrium pressure:

$$\begin{array}{cc}\hfill & {\dot{g}}_h<0\mathrm{iff}{P}_g<P_e\mathrm{and}{S}_h>0,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \mathrm{and},& {\dot{g}}_h>0\mathrm{iff}{P}_g>P_e\mathrm{and}{S}_g>0,S_w>0.\hfill \end{array}$$

(17)

At $P_g=P_e$, no reaction can occur, that is, ${\dot{g}}_h=0$, irrespective of the phase distributions. Additionally, from Equations (13)–(15), it follows that

$${\dot{g}}_h>0\Rightarrow {\dot{g}}^{C{H}_4}<0,{\dot{g}}^{H_{2}O}<0,$$

and vice-versa. In the Kim–Bishnoi kinetic model (Equation (13)), the constraints (16) and (17) are ensured through,

$$\begin{array}{ll}{k}^r>0\mathrm{and},& A_{rs}={\Gamma }_r{A}_s\\ & \mathrm{s}.\mathrm{t}.A_s=A_0{\left(1-S_h\right)}^n\hspace{1em}\mathrm{with}n>0,\\ & \begin{array}{cc}\mathrm{and}{\Gamma }_r& =\left\{\begin{array}{cc}{S}_h\hfill & \mathrm{for}\left(P_e-P_g\right)>0\hfill \\ S_g{S}_w\hfill & \mathrm{for}\left(P_e-P_g\right)\le 0,\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(18)

where $A_s$ denotes the specific surface area of the composite sediment matrix, while $A_0$ denotes the specific surface area of the primary sediment matrix. Note that in the limit of $S_h=1$ (i.e., fully clogged pores,) no reaction will occur in either direction due to unavailability of reaction surfaces. Within the scope of this work, we do not consider this limit. Hydrate dissociation is an endothermic process, and conversely, hydrate formation is an exothermic process. The heat of reaction associated with the hydrate phase change is commonly modeled as empirical functions of the form (e.g., [39]):

$$\begin{array}{c}\hfill {\dot{Q}}_h=\frac{{\dot{g}}_h}{M_h}\left(a_1+a_{2}T\right).\end{array}$$

(19)

2.2.4. Hydraulic Properties

The capillary pressure $P_c$ of the composite matrix is modeled as

$$\begin{array}{cc}\hfill & P_c=P_{c0}·f_{S_h}^{Pc}\left(S_h\right)\hfill \\ \hfill \mathrm{where} & P_{c0}=p_0{S}_{we}^{-1/\lambda }\hfill \\ \hfill \mathrm{and}, & f^{Pc}={\left(1-S_h\right)}^{-\frac{m\lambda -1}{m\lambda }}.\hfill \end{array}$$

(20)

In Equation (20), $P_{c0}$ denotes the capillary pressure of the primary matrix, and $f^{Pc}$ denotes the scaling factor which accounts for the effect of changing effective pore space due to hydrate phase change. $P_{c0}$ is parameterized using the Brooks–Corey [40] model (with soil specific parameters $p_0$ and $\lambda$, and the normalized aqueous phase saturation $S_{we}$).

The relative fluid phase permeabilities are also parameterized following the Brooks–Corey model:

$$\begin{array}{c}\hfill k_{rw}=S_{we}^{\frac{2+3\lambda }{\lambda }}\mathrm{and}{k}_{rg}={\left(1-S_{we}\right)}^2\left(1-S_{we}^{\frac{2+\lambda }{\lambda }}\right).\end{array}$$

(21)

The intrinsic permeability of the composite sediment matrix is modeled as

$$\begin{array}{cc}\hfill & K=K_0·f^K\left(S_h\right)\hfill \\ \hfill \mathrm{where}& f^K={\left(1-S_h\right)}^{\frac{5m+4}{2m}},\hfill \end{array}$$

(22)

where $K_0$ is the intrinsic permeability of the primary matrix, and $f^K$ is the scaling factor which accounts for the effect of changing effective pore space due to hydrate phase change. The scaling factors $f^{Pc}$ and $f^K$ were derived [41] based on the assumption that hydrate grows uniformly along the pore surfaces. Factor m is a measure of the sphericity of the hydrate growth. In general settings, $0<m\le 3$. For the ideal case of a spherical growth, $m=3$. The more the hydrate growth skews in the direction of the grain contacts, the lower is the m value. For example, according to the experimental investigations by [42], for hydrates formed in quartz sand, $K=K_0{\left(1-S_h\right)}^{11.4}$, implying that $m=0.225$.

2.3. Primary Variables

To solve the mathematical model numerically, we substitute the Darcy velocity (Equation (4)) and the constitutive relationships (12)–(19) in the coupled system of Equations (1)–(3) and (5). This results in a highly nonlinear system with 11 unknowns: $\left(P_g,P_w,S_g,S_w,S_h,{\chi }_w^c,{\chi }_g^{C{H}_4},{\chi }_w^{C{H}_4},{\chi }_g^{H_{2}O},{\chi }_w^{H_{2}O},T\right).$ Equation (6) provides an additional relationship for the phase saturations, reducing the number of unknown saturations to 2. Equation (7) provides a relationship for phase pressures, leaving only 1 pressure unknown. Finally, Equation (8) gives a relationship for ${\chi }_{\alpha }^{C{H}_4}$ and Equation (9) gives a relationship for ${\chi }_{\alpha }^{H_{2}O}$, thus reducing the unknown mole fractions to two. This leaves seven primary unknowns which need to be solved for the coupled system which includes four nonlinear second-order PDEs: (1), (3), and (5), one nonlinear nonhomogeneous first-order ODE (2), and two inequality constraints (11). We chose the following set of primary variables:

$$\begin{array}{c}\hfill \mathcal{P}:=\left(P_w,S_g,S_h,{\chi }_w^c,X_w^{C{H}_4},X_g^{H_{2}O},T\right).\end{array}$$

(23)

This choice of primary variables is not unique and depends on the actual application. In our case, the applications of interest arise from marine geological settings where water phase always persists but the gas phase may or may not exist (making $P_w$ a more relevant variable), and along with the hydrate phase saturations, the gas phase saturation and dissolved methane mole fraction are the most important quantities of interest, making (23) the most suitable choice for primary unknowns.

3. Numerical Solution Strategy

3.1. Space and Time Discretization of the Conservation Laws

The Equations (1)–(3) and (5) are discretized in space using a classical cell-centered finite volumes method defined on orthogonal meshes $T_h$ with $\mathcal{N}$ finite volume cells. The fluxes are evaluated using a two-point finite difference approximation of the gradients. Convective fluxes are fully upwinded. For time discretization, an implicit Euler method is used. The discretized model can be represented as a system of nonlinear algebraic equations as

$$\begin{array}{c}\hfill \mathbf{F}:=\mathbf{A}\left({\mathbf{X}}^{n+1},{\mathbf{X}}^n\right){\mathbf{X}}^{n+1}+\mathbf{B}\left({\mathbf{X}}^{n+1},{\mathbf{X}}^n\right)=\mathbf{0},\end{array}$$

(24)

where $\mathbf{X}$ denotes the solution vector which contains the discrete finite-volume approximations of the unknowns $\mathcal{P}$ at each cell center. The indices $n+1$ and n denote the solution at times $t^{n+1}$ and $t^n$.

3.2. Nonlinear Complementary Constraints

The complementarity constraints (11) can be rewritten as equivalent non-differentiable but semi-smooth functions as proposed in [15]:

$$\begin{array}{c}\hfill S_{\alpha }-max\left\{0,S_{\alpha }-\left(1-\sum _{\kappa }{\chi }_{\alpha }^{\kappa }\right)\right\}=0,\forall \alpha .\end{array}$$

(25)

These are piecewise linear with respect to the variables $S_{\alpha }$, ${\chi }_{\alpha }^{\kappa }$. Such functions are commonly referred to as complementary functions or $NCP$ functions in the literature. Some examples of other forms of such functions include the minimum function and Fischer–Burmeister function (see [29,43,44,45,46]). The complementarity constraints (11) and their equivalent form (25) are local in nature, and must hold cell-wise as

$$\begin{array}{c}\hfill \forall j\in \mathcal{N}:{\left(S_{\alpha }\right)}_j-max\left\{0,{\left(S_{\alpha }\right)}_j-\left(1-\sum _{\kappa }{\left({\chi }_{\alpha }^{\kappa }\right)}_j\right)\right\}=0,\forall \alpha .\end{array}$$

(26)

Note that the degrees of freedom of (26) can be partitioned into the following active–inactive sets:

$$\begin{array}{c}\hfill {\mathcal{A}}_{\alpha }:=\left\{j\in \mathcal{N}:{\left(S_{\alpha }\right)}_j-\left(1-\sum _{\kappa }{\left({\chi }_{\alpha }^{\kappa }\right)}_j\right)>0\right\},{\mathcal{I}}_{\alpha }:=\mathcal{N}\setminus {\mathcal{A}}_{\alpha }.\end{array}$$

(27)

The active sets ${\mathcal{A}}_{\alpha }$ correspond to the cells where phase $\alpha$ is present, while the inactive sets ${\mathcal{I}}_{\alpha }$ correspond to the cells where phase $\alpha$ is absent. Using relationships (6), (8), and (9) in Equation (26), we get the following system of non-differentiable equations:

$$\begin{array}{cc}\hfill & {\mathcal{C}}_g\left({\mathcal{P}}_j^{n+1}\right):={\left(S_g\right)}_j^{n+1}-max\left\{0,{\left(S_g\right)}_j^{n+1}-\left(1-{{\Pi }_g}_j^{n+1}{\left({\chi }_w^{C{H}_4}\right)}_j^{n+1}-{\left({\chi }_g^{H_{2}O}\right)}_j^{n+1}\right)\right\}=0\hfill \\ \hfill & {\mathcal{C}}_w\left({\mathcal{P}}_j^{n+1}\right):=1-{\left(S_g\right)}_j^{n+1}-{\left(S_h\right)}_j^{n+1}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & -max\left\{0,1-{\left(S_g\right)}_j^{n+1}-{\left(S_h\right)}_j^{n+1}-\left(1-{\left({\chi }_w^{C{H}_4}\right)}_j^{n+1}-{{\Pi }_w}_j^{n+1}{\left({\chi }_g^{H_{2}O}\right)}_j^{n+1}-{\left({\chi }_w^c\right)}_j^{n+1}\right)\right\}=0\hfill \end{array}$$

(29)

where ${\Pi }_g:={\displaystyle \frac{H_w^{C{H}_4}}{z^{C{H}_4}{P}_g}}$ and ${\Pi }_w:={\displaystyle \frac{P_g}{P_{sat}^{H_{2}O}}}$.

3.3. Semi-Smooth Newton Scheme

The system of Equations (28) and (29) is semi-smooth and piecewise differentiable. We solve these equations together with the system (24) within the same iterative loop using a generalized variant of the Newton scheme for semi-smooth problems [31]. The classical Newton method is valid in all regions where the Equations (28) and (29) are differentiable, while in other regions where Equations (28) and (29) are non-differentiable, the Jacobian can be evaluated by extending the value of the derivatives from the neighborhood of the non-differentiable regions. We approximate the Jacobian for our Newton scheme using a central difference method. To approximate the Jacobian for Equations (28) and (29), due to their piecewise smooth nature, we use the approximate active/inactive sets ${\mathcal{A}}_{\alpha }^{\left(l\right)}$ and ${\mathcal{I}}_{\alpha }^{\left(l\right)}$ at the l-th Newton step to determine the phase-wise NCP equations in each cell:

$$\begin{array}{cc}\hfill {\mathcal{C}}_g^{\left(l\right)}& =\left\{\begin{array}{cc}1-{{\Pi }_g}_j^{n+1}{\left({\chi }_w^{C{H}_4}\right)}_j^{n+1}-{\left({\chi }_g^{H_{2}O}\right)}_j^{n+1},\hfill & \hfill \mathrm{for}j\in {\mathcal{I}}_g^{\left(l\right)}\\ {\left(S_g\right)}_j^{n+1},\hfill & \hfill \mathrm{for}j\in {\mathcal{A}}_g^{\left(l\right)}\end{array}\right.\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\mathcal{C}}_w^{\left(l\right)}& =\left\{\begin{array}{cc}1-{\left({\chi }_w^{C{H}_4}\right)}_j^{n+1}-{{\Pi }_w}_j^{n+1}{\left({\chi }_g^{H_{2}O}\right)}_j^{n+1}-{\left({\chi }_w^c\right)}_j^{n+1},\hfill & \hfill \mathrm{for}j\in {\mathcal{I}}_w^{\left(l\right)}\\ 1-{\left(S_g\right)}_j^{n+1}-{\left(S_h\right)}_j^{n+1},\hfill & \hfill \mathrm{for}j\in {\mathcal{A}}_w^{\left(l\right)}.\end{array}\right.\hfill \end{array}$$

(31)

The approximate active/inactive sets ${\mathcal{A}}_{\alpha }^{\left(l\right)}$ and ${\mathcal{I}}_{\alpha }^{\left(l\right)}$ may change several times during the Newton loop, but if the Newton method converges, the final active sets will correspond to the physically correct phase state of the system. The advantage of our semi-smooth Newton scheme is that the treatment of the phase transitions is consistent within the Newton loop, which makes it easier to determine the physically correct phase state even for strongly coupled phase transitions. The Newton iteration is rather robust with respect to the initialization of the active/inactive sets, and therefore larger time step sizes can be used.

3.4. Numerical Implementation

We implemented the semi-smooth Newton scheme described in Section 3.3 for solving the nonlinear system (24), (28), and (29) within the DUNE-PDElab framework [47] which is based on C++. For solution of the linear system arising from the Newton linearization, we used the SUPERLU linear solver [48]. For parallel computations, we used a parallel algebraic multigrid (AMG) solver which uses a stabilized bi-conjugate gradient method as a preconditioner and a symmetric successive over-relaxation smoothening algorithm. The AMG solver is built into the dune-istl library (https://www.dune-project.org/modules/dune-istl/). For making the numerical computations, we used the NEC HPC-Linux-Cluster which is part of a hybrid NEC high-performance system at the University Computing Centre of the Christian Albrecht Universität, Kiel, Germany.

4. Numerical Examples

Here, we present two numerical examples based on the geological setting of the Black Sea. In Example 1, we simulated the sedimentation-driven gas migration through the GHSZ in the Black Sea in a 1D domain. We compared the performance of our semi-smooth Newton scheme against a PVS scheme to show the robustness of our numerical scheme. In Example 2, we simulated a gas production scenario in a 2D domain. The hydrate phase was randomly distributed within the hydrate layer with saturations ranging between 0.0 and 0.6. Through this example we show that our numerical scheme can robustly handle phase transitions even when the phase distributions and the corresponding permeability and porosity profiles are highly complex with large variations.

4.1. Example 1: Gas Migration through Gas Hydrate Stability Zone (GHSZ)

In this example, we simulated the gas hydrate dynamics driven by the changes in temperature, pressure, and salinity conditions as a result of the sediment deposition in the highly dynamic geological setting of the Black Sea. Through this field-scale environmental application, we aimed to demonstrate the complexities and challenges associated with the highly coupled phase transitions in natural gas hydrate systems, and show the robustness of our semi-smooth Newton scheme in realistic settings. We also compared the performance of our semi-smooth Newton scheme with that of a primary variable switching scheme.

4.1.1. Problem Setting

The geological setting for this problem was based on the Danube paleo delta which consists of stacked channel-levee systems that were active during glacial times when the water level was approximately 100–150 m lower than today [49]. For our problem, we chose a buried channel-levee (BCL) complex (blue color, Figure 2) west of the Viteaz canyon, the main Danube paleo channel, which has buried the BCL over the past 75 ka (green color, Figure 2). The BCL is believed to have deposited its levees essentially in two main events correlating to oxygen isotope stages 8 and 6 [49], that is, between 320 ka and 75 ka BP (brown color, Figure 2 [50]). These two active phases of the BCL correspond to limnic stages of the Black Sea that have been documented, for example by low sulfur contents, in the sedimentary record of DSDP drill Site 379 in the eastern basin of the Black Sea [51]. For the past 500 ka, this drill core identifies five marine stages that are interrupted by four limnic stages—that is, intervals when the sea level dropped below the depth of Bosphorus sill (today at 40 m water depth), thereby separating the Black Sea from saltwater inflow from the Sea of Marmara and the Mediterranean Sea.

In our problem, we were interested in simulating how the deposition of the brown and green sediment layers affects the gas hydrate stability zone (GHSZ) that was established 300 ka BP, that is, in the blue sediments (Figure 2). Hence, we based the initial setting $\left(\mathrm{at}\mathrm{time}t=t_0=0\right)$ on the paleo conditions existing at 300 ka BP. We chose an arbitrary $1D$ segment $A-B$ located in the eastern levee as our computational domain (see Figure 2). Point A of the computational domain corresponds to the paleo seafloor at 300 ka BP (PSF-C) (i.e., $z=z_A=0$ m, where z denotes the depth below the sea floor). Point B is located at $z=z_B=-800$ m. At $t=t_0$, we assumed a hydrostatic pressure at point A of ${\left.P_w\right|}_{z=0,t=0}=P_{sf}=15$ MPa, corresponding to a water depth of roughly 1500 m, and a bottom water temperature of ${\left.T\right|}_{z=0,t=0}=T_{sf}=4$${}^{\circ }$C, corresponding to glacial conditions in the Black Sea [50]. We assumed that the initial pressure distribution within the computational domain follows a hydrostatic gradient, and the initial temperature distribution follows a steady-state geothermal gradient of 35 ${}^{\circ }\mathrm{C}/\mathrm{km}$. The initial conditions for all the primary variables are listed in

Table 1. Initial and bounday conditions for Example 1 (Section 4.1).

Initial Conditions
at $t=0$, and 0 m $\ge z\ge$ $-800$ m	$P_w$	=	$P_{sf}+{\rho }_{w}g\left(z_{sf}-z\right)$
	where $z_{sf}=0$ m is the sea floor,
	and $P_{sf}=15$ MPa is the water pressure at the sea floor.
	T	=	$T_{sf}+d_z{T}_G\left(z_{sf}-z\right)$
	where $T_{sf}=4$ ${}^{\circ }$C is the bottom water temperature,
	and $d_z{T}_G=35$ ${}^{\circ }\mathrm{C}/\mathrm{km}$ denotes the regional geothermal temperature gradient.
	$S_g$	=	0
	${\chi }_w^{C{H}_4}$	=	0
	${\chi }_g^{H_{2}O}$	=	${\chi }_{g,sat}^{H_{2}O}\left({\left.P_g\right|}_{t=0},{\left.T\right|}_{t=0}\right)$
	${\chi }_w^c$	=	$5.5\mathrm{mmol}/\mathrm{mol}$
at $t=0$,
if,    $-320$ m $\ge z\ge$ $-400$ m,	$S_h$	=	$0.3\left({\displaystyle \frac{400+z}{400-320}}\right)\left({\displaystyle \frac{z+320}{400-320}}\right)$
else if,    $z>-320$ m or $z<-400$ m	$S_h$	=	0
Boundary Conditions
$t>0$, and $z=0$ m	$P_w$	=	$P_{sf}+{\rho }_{s}g{v}_{s,z}\left(t_n+\Delta t\right)$
	T	=	$T_{sf}+d_z{T}_G{v}_{s,z}\left(t_n+\Delta t\right)$
	$S_g$	=	0
	${\chi }_w^c$	=	${\left.{\chi }_w^c\right|}_{t=0}$
$t>0$, and $z=-800$ m	${\partial }_z{P}_w$	=	0
	${\partial }_{z}T$	=	$d_z{T}_G$
	$v_{g,z}$	=	0
	${\partial }_z{\chi }_w^c$	=	0

. Based on the initial pressure, temperature, and salinity conditions, we could estimate the location of the base of the GHSZ (bGHSZ), that is, the point of intersection of the gas phase pressure and the gas hydrate equilibrium pressure curves plotted along the sediment depth. The gas hydrates are stable above bGHSZ where $P_g\ge P_{eq}$, and unstable below. (see Figure 3a). For this setting, the initial bGHSZ was located 400 m below point A, and we assumed that a hydrate layer of 80 m thickness and 30% peak saturation was located just above this initial bGHSZ. For the sake of simplicity, we assumed that the deposition of the brown and green sediment layers occurred continuously over 300,000 years at a constant sedimentation rate of $v_{s,z}=0.1$ cm/year. This does not reflect the true depositional history, but rather simulates a scenario of a low average sedimentation rate. Figure 3b shows schematically how the sedimentation shifts the GHSZ. Basically, due to the sedimentation process, the sea floor rises over time. At any time $t=t_n>t_0$, the corresponding sea floor PSF-n is located at $z=z_{\mathrm{PSF}-\mathrm{n}}=z_A+v_{s,z}{t}_n$, and within a time increment $\Delta t$, a new sediment layer of thickness $\Delta z=v_{s,z}\Delta t$ is deposited on top of PSF-n. We assumed that $\Delta t$ was small enough for temperature and pressure to reach a steady state within the new sediment layer. The pressure and temperature at any sea floor PSF-n were are fixed at ${\left.P_w\right|}_{z_{\mathrm{PSF}-\mathrm{n}},t_n}=P_{sf}$ and ${\left.T\right|}_{z_{\mathrm{PSF}-\mathrm{n}},t_n}=T_{sf}$, respectively. Note that here we ignored any changes in sea level and bottom water temperature during the geological past. Due to the sedimentation, the temperature and pressure at the top boundary of our computational domain (i.e., point A at $z=0$ m) increased over time, which in turn shifted the base of the GHSZ upwards. (refer to Table 1 for a list of the boundary conditions, and

Table 2. Material properties and model parameters for Examples 1 and 2.

Property ${}^{1,2}$			Example 1	Example 2
Water
density	${\rho }_w$	kg/m${}^3$	$1030.21$	$1027+0.45P_w\left[\mathrm{MPa}\right]-0.15\left(T\left[{}^{\circ }\mathrm{C}\right]-10\right)+0.3521\times {10}^3\left({\chi }_w^c-0.0096\right)$
dynamic viscosity	${\mu }_w$	Pa.s	$0.00136$	$0.001792exp\left(\phantom{{\left({\displaystyle \frac{}{}}\right)}^0}-1.94\right.-4.80\left({\displaystyle \frac{273.15}{T}}\right)+\left.6.75{\left({\displaystyle \frac{273.15}{T}}\right)}^2\right)$
thermal conductivity	$k_w^{th}$	W/m/K	$0.59$	$0.57153\left(1\right.+0.003T\left[{}^{\circ }\mathrm{C}\right]-1.025\times {10}^{-5}{\left(T\left[{}^{\circ }\mathrm{C}\right]\right)}^2+6.53\times {10}^{-10}{P}_w-\left.0.0797{\chi }_w^c\right)$
specific heat capacity	$C{p}_w$	J/kg/K	3945
saturation vapor pressure	$P_{sat}^{H_{2}O}$	Pa	$Pcexp\left({\displaystyle \frac{1}{Tr}}\left(c_1(1-Tr)\right.\right.+c_2{(1-Tr)}^{1.5}+c_3{(1-Tr)}^3+c_4{(1-Tr)}^{3.5}+c_5{(1-Tr)}^4+\left.\left.c_6{(1-Tr)}^{7.5}\right)\phantom{{\displaystyle \frac{}{}}}\right)$
			where $Pc=22.064$ MPa, $Tc=647.096$ K, $Tr=T/Tc$, and $c_1=-7.85951783$, $c_2=1.84408259$, $c_3=-11.7866497$, $c_4=22.6807411$, $c_5=-15.9618719$, $c_6=1.80122502$.
diffusion coefficient	$D_g^{H_{2}O}$	m${}^2$/s	$0.638\times {10}^{-6}$	$2.26\times {10}^{-9}T+{\displaystyle \frac{0.002554}{P_g}}$
Methane
density	${\rho }_g$	kg/m${}^3$	$0.002756{\displaystyle \frac{P_g}{T}}$	${\displaystyle \frac{P_g}{z_{C{H}_4}{R}_{C{H}_4}T}}$    where $R^{C{H}_4}=\frac{8314.5}{16.04}$, and $z^{C{H}_4}$ is estimated using Peng Robinson EoS [52].
dynamic viscosity	${\mu }_g$	Pa.s	$1.4055\times {10}^{-5}$	${\mu }_0\left({\displaystyle \frac{273.15+162}{T+162}}\right){\left({\displaystyle \frac{T}{273.15}}\right)}^{1.5}$
				where ${\mu }_0=1.0707\times {10}^{-5}-4.8134\times {10}^{-14}{P}_g-4.1719\times {10}^{-20}{P}_g^2+7.3232\times {10}^{-28}{P}_g^3$
thermal conductivity	$k_g^{th}$	W/m/K	$0.03121$	$a_0+a_{1}T+a_2{T}^2+a_3{T}^3$
				where $a_0=-0.008863$, $a_1=0.000242$, $a_2=-0.6997\times {10}^{-6}$, and $a_3=0.1225\times {10}^{-8}$
specific heat capacity	$C{p}_g$	J/kg/K	$2168.65$	$1238+3.13T+7.905\times {10}^{-4}{T}^2-6.858\times {10}^{-7}{T}^3$
solubility constant	$H_w^{C{H}_4}$	Pa	$exp\left(log\left(P_{sat}^{H_{2}O}\right)-{\displaystyle \frac{11.0094}{Tr}}+4.8362{\displaystyle \frac{{(1-Tr)}^{0.355}}{Tr}}+12.5220exp(1-Tr)T{r}^{-0.41}\right)$
diffusion coefficient	$D_w^{C{H}_4}$	m${}^2$/s	$1.57\times {10}^{-11}$	$1.57\times {10}^{-11}\left({\displaystyle \frac{P_w}{1.0135\times {10}^5}}\right)exp\left(-{\displaystyle \frac{0.003475}{T}}\right)$
Hydrate
density	${\rho }_h$	kg/m${}^3$	920
hydration number	$N_h$	−	$5.90$
thermal conductivity	$k_h^{th}$	W/m/K	$0.5$
specific heat capacity	$C{p}_h$	J/kg/K	2327	$\left(1.937T^3\right.-1.5151T^2+3.9554T-\left.342.7\right)\times {10}^3$
Salt
diffusion coefficient	$D_w^c$	m${}^2$/s	${10}^{-9}$
Soil
density	${\rho }_s$	kg/m${}^3$	2600
thermal conductivity	$k_s^{th}$	W/m/K	$3.0$
specific heat capacity	$C{p}_s$	J/kg/K	1000
Hydrate Phase Change Kinetics
hydrate equilibrium pressure	$P_e$	Pa	$exp\left(38.592-{\displaystyle \frac{8533.8}{T}}+4.4824{\chi }_w^c\right)$
kinetic rate constant	$k^r$	mol/m${}^2$/Pa/s	${10}^{-17}$	${10}^{-12}$
specific surface area	$A_0$	m${}^2$/m${}^3$	${10}^5$
heat of reaction	${\dot{Q}}_h$	W/m${}^3$	${\displaystyle \frac{{\dot{g}}_h}{M_h}}\left(56599-16.744T\right)$
Hydraulic Properties
absolute intrinsic permeability	$K_0$	m${}^2$	${10}^{-15}$	${10}^{-13}$
total porosity	$\varphi$	−	$0.5$	$0.3$
Brooks–Corey parameters	$p_0$, $\lambda$	Pa,−	$5\times {10}^4$, $1.2$
sphericity parameter	m	−	1
residual saturations	$S_{wr}$, $S_{gr}$	−,−	0,0

¹ See [53] for references to $P_{sat}^{H_{2}O}$ an $H_w^{C{H}_4}$; [20] for references to $k_w^{th}$, $C{p}_w$, ${\mu }_{0,g}$, ${\rho }_h$, $N_h$, $k_h^{th}$, $C{p}_h$, $D_w^c$, ${\rho }_s$, $k_s^{th}$, and $C{p}_s$; and [54] for references to ${\mu }_w$, $D_g^{H_{2}O}$, ${\mu }_g$, $D_w^{C{H}_4}$, $k_g^{th}$, $C{p}_g$, and ${\dot{Q}}_h$. ² To evaluate the property values for Example 1, we considered a reference state of $T=9$ ${}^{\circ }$C, $P=10\mathrm{MPa}$, and ${\chi }_w^c=5.5$ mmmol/mol.

for a list of material properties and parameters).

The main challenge in simulating this setting is that, as the hydrate layer gets buried below the GHSZ due to ongoing sedimentation, it starts to dissociate from the bottom, and the gas phase appears in a narrow region below the GHSZ. The saturation of the free gas phase is very small, typically less than 5%. The gas migrates upwards due to its buoyancy, but since the overlying hydrate layer has a much lower permeability, the free gas tends to pool below the region where the hydrate saturation is the highest, thereby building up the pore pressure. The local dilution of the pore water salinity due to fresh water release and the local cooling effect due to hydrate dissociation also give strong feedbacks to both the hydrate equilibrium pressure, as well as the solubility of the gas in the aqueous phase. These competing effects often cause the mathematical model to rapidly switch back and forth between a single phase model and a two-phase model with respect to the $C{H}_4-H_{2}O$ system, especially when the gas phase appears in the domain for the first time.

4.1.2. Numerical Simulation and Results

The computational domain was discretized uniformly into 1600 finite volumes along the Z-axis. The maximum time step size was chosen as $d{t}_{max}=10$ years. An adaptive time-stepping strategy was used where the time step size was controlled heuristically based on the number of Newton iterations per time integration step. If the number of Newton iterations was more than ${\ell }_h$, the time step size for the next time integration step was decreased by 25%, whereas if the number of Newton iterations was less than ${\ell }_l$, the time step size was increased by 10%. Between ${\ell }_l$ and ${\ell }_h$ Newton iterations, $dt$ was not changed. The choice of ${\ell }_l$ and ${\ell }_h$ depends on the problem setting and, as a general rule, in our numerical scheme we considered ${\ell }_l\ge 4$ and ${\ell }_h={\ell }_l+4$. In this example, we chose ${\ell }_l=8$ and ${\ell }_h=12$.

We performed the numerical simulations with our semi-smooth Newton (NCP) scheme, and for comparison, also with the more common primary variable switching (PVS) approach of [11]. We implemented this PVS scheme within the same software framework as our NCP approach (i.e., DUNE-PDElab, version 2.6.0; https://www.dune-project.org/modules/dune-pdelab/). For both schemes, we used the same discretization scheme (Section 3.1) and the same linear solver (SuperLU). For the Newton solver, we used the same convergence criteria, and for the adaptive time-stepping strategy, we used identical control parameters. Note that we implemented only a sequential version of the PVS scheme. Therefore, for those examples where we compared the solution of our NCP scheme with the PVS scheme, we performed all numerical simulations only in a sequential mode.

In Figure 4a, we can see that the PVS scheme needed ≈50% more CPU time to solve only until ≈45% of the problem time even though the PVS scheme needed less time per calculation due to fewer degrees of freedom compared to the NCP scheme.

In Figure 5, the snapshots of $S_h$, $S_g$, ${\chi }_w^{C{H}_4}$, and ${\chi }_w^c$ are plotted at three times: $t_1$ = 22,500 years, $t_2$ = 135,000 years, and $t_3$ = 300,000 years. Time $t_1$ corresponds to the instant when the gas phase first appears in the domain. We can see that both the PVS and the NCP schemes are in agreement about where the gas phase appears and in what saturation. Time $t_2$ corresponds to the time up to which the PVS scheme could solve in 240 CPU-hours (at which point the simulation was aborted due to the large run time). The results of the PVS and the NCP schemes match very well, and serve as a good validation for our numerical implementation. Time $t_3$ corresponds to the end-time for this problem. Only the solution of the NCP scheme is plotted for this time step. The base of the GHSZ shifts upwards by 300 m due to sedimentation over a period of 300,000 years. The results show that as the GHSZ risis towards the sea floor, the hydrate layer dissociates, generating methane below the base of the GHSZ. Through a combination of dissolution, diffusion, and buoyancy effects, methane transports into the GHSZ where the gas hydrate layer re-forms. The hydrate layer follows the base of the GHSZ, but shrinks along the way as more and more gas dissolves and diffuses away. In Figure 6, this process of hydrate dissociation → gas migration → hydrate reformation is shown in greater detail by zooming in on the time axis between 60,000 years $\le t\le$ 120,000 years. These processes are quite complex and nonlinear. As the gas hydrate dissociates from the bottom, the free gas rises upwards with a decreasing speed due to a decreasing permeability in the hydrate phase. When the gas phase crosses the region with maximum $S_h$, the speed of gas migration starts to increase until it escapes the hydrate layer on the other side, where a new hydrate layer starts to form, as shown in Figure 6c. This new hydrate layer continues to grow using the free gas supplied by the dissociation of the old hydrate layer, as shown in Figure 6c–e.

In Figure 4b, the evolution of $dt$ is plotted over the problem time for both schemes. We can see that at $t=22,500$ years, when the gas phase first appears, $dt$ breaks down for both schemes. However, the reduction of $dt$ for NCP scheme was not as severe as that for the PVS scheme. The time step size gradually recovers as the gas slowly migrates upwards through the hydrate layer, but breaks down again around $t=60,000$ years, when the free gas crosses the region with peak $S_h$.

4.2. Example 2: Gas Production through Depressurization

In this example, we simulated a gas production scenario where the gas hydrate reservoir was destabilized through depressurization. We considered a single-well configuration, and based the model parameters, material properties, and initial conditions within the reservoir on the geological setting of the Black Sea. The objective of this example was to demonstrate the robustness of our scheme in handling complex phase transitions even in those settings where the reaction rates are large, the hydrate phase distributions are highly heterogeneous, and the permeability typically varies over two-to-four orders of magnitude. Such settings are commonly found in natural gas hydrate systems which occur in turbidite formations containing thin hydrate layers sandwiched between thin layers of silty to clayey sediments.

4.2.1. Problem Setting

We considered a $2D$ axisymmetric domain, $\Omega$, having dimensions 1000 m × 1000 m, as shown in Figure 7. The sea floor, $\partial {\Omega }_{sf}$, was prescribed at $z=0$ m. The depressurization well, $\partial {\Omega }_{well}$, was located at $r=0$ m, 0 m $\ge z\ge$$-400$ m. The bottom water temperature at the sea floor was $T_{sf}=4$ ${}^{\circ }$C, and the hydrostatic pressure at the sea floor was $P_{sf}=15$ MPa. We assumed that the absolute intrinsic permeability and the total porosity of the primary soil skeleton were $K={10}^{-13}{\mathrm{m}}^2$ and $\varphi =0.3$, respectively. At $t=0$, we assumed that the domain was fully saturated with saline water, and there was no free gas phase in the domain. Additionally, the aqueous phase contained no dissolved methane. For the aqueous phase pressure, we considered a hydrostatic pressure gradient, and for the temperature, we prescribed a regional geothermal gradient along the depth, $d_z{T}_G=35$${}^{\circ }$C/km. The bGHSZ was located at $z=-400$ m. We considered that an 80-m-thick gas hydrate layer, ${\Omega }_H$, existed right above the bGHSZ. To show the robustness of our scheme with respect to complex phase transitions, we prescribed a random distribution of the hydrate phase within this layer, s.t. the hydrate saturation ranged from $0.0$ to $0.6$, and the corresponding absolute permeability ranged from ${10}^{-13}{\mathrm{m}}^2$ to $1.6\times {10}^{-15}{\mathrm{m}}^2$. For $t>0$, a pressure of ${\left.P_w\right|}_{\partial {\Omega }_{well}}=8$ MPa was prescribed at the production well to simulate gas production through depressurization. At the sea floor, the temperature, pressure, and salinity conditions remained constant and equal to the initial values. At the bottom boundary, $\partial {\Omega }_B$, the regional geothermal gradient was maintained. The initial and boundary conditions are listed in

Table 3. Initial and bounday conditions for Example 2 (Section 4.2).

Initial Conditions
tag: $\Omega$ at $t=0$, 0 m $\le r\le$ 1000 m and 0 m $\ge z\ge$ $-1000$ m	$P_w$	=	$P_{sf}+{\rho }_{w}g\left(z_{sf}-z\right)$
	where $z_{sf}$ denotes the sea floor, $z_{sf}=0$ m,
	and $P_{sf}$ denotes the water pressure at the sea floor, $P_{sf}=15$ MPa.
	T	=	$T_{sf}+d_z{T}_G\left(z_{sf}-z\right)$
	where $T_{sf}=4$ ${}^{\circ }$C denotes the temperature at the sea floor,
	and, $d_z{T}_G=35{}^{\circ }\mathrm{C}/\mathrm{km}$ denotes the regional geothermal temperature gradient.
	$S_g$	=	0
	${\chi }_w^{C{H}_4}$	=	0
	${\chi }_g^{H_{2}O}$	=	${\chi }_{g,sat}^{H_{2}O}\left({\left.P_g\right|}_{t=0},{\left.T\right|}_{t=0}\right)$
	${\chi }_w^c$	=	$5.5\mathrm{mmol}/\mathrm{mol}$
at $t=0$, and  0 m $\le r\le$ 1000 m,
tag: ${\Omega }_H$:       $-320$ m $\ge z\ge$ $-400$ m	$S_h$	=	rand $\left[0,0.6\right]$
tag: $\Omega -{\Omega }_H$:    $z>-320$ m or $z<-400$ m	$S_h$	=	0
Boundary Conditions
tag: $\partial {\Omega }_{well}$	$P_w$	=	$8\mathrm{MPa}$
$t>0$,    $r=0$ and 0 m $\ge z\ge$ $-400$ m	$\nabla T$	=	0
	${\mathbf{v}}_g$	=	0
	$\nabla {\chi }_w^c$	=	0
tag: $\partial {\Omega }_{sf}$	$P_w$	=	$P_{sf}$
$t>0$,    $z=0$ and 0 m $\le r\le$ 1000 m	T	=	$T_{sf}$
	$S_g$	=	0
	${\chi }_w^c$	=	${\left.{\chi }_w^c\right|}_{t=0}$
tag: $\partial {\Omega }_R$	${\mathbf{v}}_w$	=	0
$t>0$,    $r=1000$ m and 0 m $\ge z\ge$ $-1000$ m	$\nabla T$	=	0
	${\mathbf{v}}_g$	=	0
	$\nabla {\chi }_w^{c}g$	=	0
tag: $\partial {\Omega }_B$	${\mathbf{v}}_w$	=	0
$t>0$,    $z=-1000$ m and 0 m $\le r\le$ 1000 m	$\nabla T$	=	$\nabla T_G$
	${\mathbf{v}}_g$	=	0
	$\nabla {\chi }_w^{c}g$	=	0
tag: $\partial {\Omega }_L$	${\mathbf{v}}_w$	=	0
$t>0$,    $r=0$ m and $-400$ m $>z\ge$ $-1000$ m	$\nabla T$	=	0
	${\mathbf{v}}_g$	=	0
	$\nabla {\chi }_w^{c}g$	=	0

, and the material properties and model parameters are listed in Table 2.

4.2.2. Numerical Simulation and Results

The computational domain was discretized into a total of 20,276 quadrilateral elements. The gas production process was simulated until $t_{end}=360$ days. The maximum time step size was chosen as $d{t}_{max}$ = 36,000 s, and the time step size $dt$ was controlled adaptively using the heuristic strategy discussed in Section 4.1.2 with ${\ell }_l=4$ and ${\ell }_h=8$. The simulation was run in parallel on four processing units and required a total of 20 CPU-hours. We identified a domain of interest, ${\Omega }_I$, for the gas production process as: 0 m $\le r\le$ 250 m, $-150$ m $\ge z\ge$$-550$ m. Outside this domain, pressure, temperature, and saturation profiles did not change much. However, the large size of the domain is necessary to ensure that effects of depressurization do not reach $\partial {\Omega }_R$, and the geothermal gradient is maintained at $\partial {\Omega }_B$.

The main quantities of interest (QoIs) for gas production in gas hydrate reservoirs are $S_h$, $S_g$, ${\chi }_w^{C{H}_4}$, and the GHSZ. The snapshots of the QoIs within ${\Omega }_I$ are plotted in Figure 8 for times $t=10$ days, $t=90$ days, and $t=360$ days. We can see that over a period of one year, roughly 100 m of the reservoir was effectively depressurized, and the hydrate layer was fully dissociated within a zone of roughly 15 m around the production well. Due to the relatively high pore pressures, most of the methane was produced from the aqueous phase, and the saturation of the free gas phase remained well below 10%.

We compared the performance of the numerical scheme for this example with that of a reference gas production test case. The setting of the reference test case was the same as described in Section 4.2.1, except that the hydrate distribution in ${\Omega }_H$ was homogeneous and had a uniform saturation of $0.6$. A snapshot of the QoIs in ${\Omega }_H$ for the reference test case is plotted in Figure 9 at time step $t=360$ days. By comparing the behavior of the numerical solution with the reference test case, we can ensure that the random hydrate distribution did not introduce any artificial numerical artifacts in the numerical solution. In Figure 10, we can also see that despite the random distribution of the hydrate phase and the large variations in the sediment permeability, the semi-smooth Newton scheme was able to handle the phase transitions quite robustly without significant loss in performance, as indicated by the evolution of the time-step sizes.

5. Conclusions

In this article, we presented a mathematical model for non-isothermal multi-phase multi-component reactive transport processes in methane hydrate reservoirs. The methane hydrate phase transitions were modeled as a non-equilibrium based kinetic process, and the phase transitions within the ${CH}_4-H_{2}O$ system were modeled as a VLE process. The inequality conditions resulting from the VLE assumption were cast as Kharush–Kuhn–Tucker (KKT) equality conditions which were implemented within a semi-smooth Newton scheme using an active-set strategy. Note that in the context of gas hydrate models, a similar nonlinear complementary constraints approach was also used by [55] to develop a semi-smooth Newton strategy for solving a Stefan-type problem involving equilibrium-based hydrate phase transition.

In many widely used multi-phase multi-component gas hydrate reservoir simulators, PVS is the method of choice for handling phase transitions and the vanishing gas phase. In general, the PVS method has the advantage that the numerical model has fewer degrees of freedom, and therefore can perform numerical calculations faster. However, in the case of gas hydrate models, this advantage is most often lost because the phase transitions occur at different rates, and are highly coupled and extremely sensitive to the changes in the local thermodynamic state, which leads to rapidly switching phase states. The PVS scheme often shows poor convergence and requires much smaller time steps. We demonstrated this in Example 1 (Section 4.1), where we simulated gas migration through the GHSZ in the highly dynamic geological setting of the Black Sea over paleo time-scales. Note in Table 2 that, for Example 1, we greatly simplified the problem setting by neglecting the functional dependence of the material properties (e.g., densities, viscosities, heat capacities, and thermal conductivities) on local temperature, pressure, and salinity conditions, thereby reducing the nonlinearities, and we considered only a $1D$ setting with homogeneous phase distributions across the domain. Despite these simplifications, the PVS scheme performed much more poorly compared to our semi-smooth Newton scheme. In Example 2 (Section 4.2), we considered another very important application of methane hydrate models, viz., gas production through depressurization. We considered a random distribution of the hydrate phase (and consequently, a random permeability field and capillary pressures ranging over multiple orders of magnitude). We also included strongly nonlinear functional dependencies of the material properties on the local thermodynamic state (see Table 2). Through this example we demonstrated that our semi-smooth Newton scheme can robustly handle even very complex field-scale problems.