An all-at-once Newton strategy for methane hydrate reservoir models

TL;DR

This paper introduces a novel all-at-once Newton method with nonlinear complementary constraints for modeling methane hydrate reservoirs, improving convergence and robustness over traditional phase transition handling methods.

Contribution

The authors develop a nonlinear complementary constraints approach integrated with a non-smooth Newton scheme to better handle phase transitions in hydrate reservoir models.

Findings

Enhanced convergence in highly nonlinear phase transition problems

Robustness demonstrated on field-scale Black Sea applications

Outperforms traditional primary variable switching methods

Abstract

Marine gas hydrate systems are characterized by highly dynamic transport-reaction processes in an essentially water-saturated porous medium that are coupled to thermodynamic phase transitions between solid gas hydrates, free gas and dissolved methane in the aqueous phase. These phase transitions are highly nonlinear and strongly coupled, and cause the mathematical model to rapidly switch the phase states and pose serious convergence issues for the classical Newton's method. One of the common methods of dealing with such phase transitions is the primary variable switching (PVS) method where the choice of the primary variables is adapted locally to the phase state **outside** the Newton loop. In order to ensure that the phase states are determined accurately, the PVS strategy requires an additional iterative loop, which can get quite expensive for highly nonlinear problems. For methane hydrate reservoir models, the PVS method shows poor convergence behaviour and often leads to extremely small time step sizes. In order to overcome this issue, we have developed a nonlinear complementary constraints method (NCP) for handling phase transitions, and implemented it within a non-smooth Newton's linearization scheme using an active-set strategy. Here, we present our numerical scheme and show its robustness through field scale applications based on the highly dynamic geological setting of the Black Sea.

Tables

Table 1. Table 1: Initial and bounday conditions for the Example 1 (Sec. 4.1 ).

Initial conditions
at $t=0$ , and $𝐱\in \Omega$	$P_w$	$=$	$2\text{ MPa}$
	$S_g$	$=$	$0$
	$S_h$	$=$	$0.3$
	${\chi }_w^{C{H}_4}$	$=$	$0$
	${\chi }_g^{H_{2}O}$	$=$	${\chi }_{g,sat}^{H_{2}O}({P_g|}_{t=0},{T|}_{t=0})$
	${\chi }_w^c$	$=$	$5.5\text{ mmol/mol}$
	$T$	$=$	$4^{o}C$
Boundary conditions
for $𝐱\in \partial \Omega$,	$P_w$	$=$	$2\text{ MPa}$
	${𝐯}_w$	$=$	$0$
	$P_w$	$=$	$5\text{ MPa}$
$t\ge 450\text{ hrs}$	$P_w$	$=$	$5\text{ MPa}+10.\times (t-450\times 3600)$
for $𝐱\in \partial \Omega$, and $t>0$	${𝐯}_g$	$=$	$0$
	$\nabla {\chi }_w^c$	$=$	$0$
	$T$	$=$	$4^{o}C$

Table 2. Table 2: Initial and bounday conditions for Example 2 (Sec. 4.2 ).

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^0$C is the bottom water temperature,
	and, $d_z{T}_G={35}^0\text{C/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}({P_g|}_{t=0},{T|}_{t=0})$
	${\chi }_w^c$	$=$	$5.5\text{ mmol/mol}$
at $t=0$,
if,  $-320$m $\ge z\ge$ $-400$m ,	$S_h$	$=$	$0.3\left(\frac{400+z}{400-320}\right)\left(\frac{z+320}{400-320}\right)$
else if,  $z>-320$m or 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$	$=$	${{\chi }_w^c|}_{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$

Table 3. Table 3: Initial and bounday conditions for Example 3 (Sec. 4.3 ).

Initial conditions
tag: $\Omega$	$P_w$	$=$	$P_{sf}+{\rho }_{w}g\left(z_{sf}-z\right)$
at $t=0$ ,	where, $z_{sf}$ denotes the sea floor, $z_{sf}=0$m,
$0$m $\le r\le$ $1000$m and $0$m $\ge z\ge$ $-1000$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^0$C denotes the temperature at the sea floor,
	and, $d_z{T}_G={35}^0\text{C/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}({P_g|}_{t=0},{T|}_{t=0})$
	${\chi }_w^c$	$=$	$5.5\text{ mmol/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$[0,0.6]$
tag: $\Omega -{\Omega }_H$:  $z>-320$m or m	$S_h$	$=$	$0$
Boundary conditions
tag: $\partial {\Omega }_{well}$	$P_w$	$=$	$8\text{ MPa}$
$t>0$,  $r=0$ and $0$m $\ge z\ge$ $-400$m	$\nabla T$	$=$	$0$
	${𝐯}_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$	$=$	${{\chi }_w^c|}_{t=0}$
tag: $\partial {\Omega }_R$	${𝐯}_w$	$=$	$0$
$t>0$,  $r=1000$m and $0$m $\ge z\ge$ $-1000$m	$\nabla T$	$=$	$0$
	${𝐯}_g$	$=$	$0$
	$\nabla {\chi }_w^{c}g$	$=$	$0$
tag: $\partial {\Omega }_B$	${𝐯}_w$	$=$	$0$
$t>0$,  $z=-1000$m and $0$m $\le r\le$ $1000$m	$\nabla T$	$=$	$\nabla T_G$
	${𝐯}_g$	$=$	$0$
	$\nabla {\chi }_w^{c}g$	$=$	$0$
tag: $\partial {\Omega }_L$	${𝐯}_w$	$=$	$0$
$t>0$,  $r=0$m and $-400$m $>z\ge$ $-1000$m	$\nabla T$	$=$	$0$
	${𝐯}_g$	$=$	$0$
	$\nabla {\chi }_w^{c}g$	$=$	$0$

Table 4. Table 4: Material properties and model parameters for Examples 1, 2, and 3.

Property${}^{\text{2 3}}$			Example 1	Example 2	Example 3
water
density	${\rho }_w$	$kg/m^3$	$1026.77$	$1030.21$	$1027+0.45P_w\left[MPa\right]-0.15\left(T\left[{}_{}{}^{0}C\right]-10\right)+0.3521\times {10}^3\left({\chi }_w^c-0.0096\right)$
dynamic viscosity	${\mu }_w$	$Pa.s$	$0.0014$	$0.00136$	$0.001792\exp \left(-1.94-4.80\left(\frac{273.15}{T}\right)+6.75{\left(\frac{273.15}{T}\right)}^2\right)$
thermal conductivity	$k_w^{th}$	$W/m/K$	$0.5854$	$0.59$	$0.57153\left(1+0.003T\left[{}_{}{}^{o}C\right]-1.025\times {10}^{-5}{\left(T\left[{}_{}{}^{o}C\right]\right)}^2+6.53\times {10}^{-10}{P}_w-0.0797{\chi }_w^c\right)$
specific heat capacity	$C{p}_w$	$J/kg/K$	$3945$
saturation vapour pressure	$P_{sat}^{H_{2}O}$	$Pa$	$1072.92$	$Pc\exp \left(\frac{1}{Tr}\left(c_1\left(1-Tr\right)+c_2{\left(1-Tr\right)}^{1.5}+c_3{\left(1-Tr\right)}^3+c_4{\left(1-Tr\right)}^{3.5}+c_5{\left(1-Tr\right)}^4+c_6{\left(1-Tr\right)}^{7.5}\right)\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.637\times {10}^{-6}$	$0.638\times {10}^{-6}$	$2.26\times {10}^{-9}T+\frac{0.002554}{P_g}$
methane
density	${\rho }_g$	$kg/m^3$	$19.605$	$0.002756\frac{P_g}{T}$	$\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 PengRobinson1970 .
dynamic viscosity	${\mu }_g$	$Pa.s$	$1.1045\times {10}^{-5}$	$1.4055\times {10}^{-5}$	${\mu }_0\left(\frac{273.15+162}{T+162}\right){\left(\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.03107$	$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$	$2165.24$	$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$	$1,343\times {10}^{11}$	$\exp \left(\log \left(P_{sat}^{H_{2}O}\right)-\frac{11.0094}{Tr}+4.8362\frac{{\left(1-Tr\right)}^{0.355}}{Tr}+12.5220\exp \left(1-Tr\right)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}$	$1.57\times {10}^{-11}\left(\frac{P_w}{1.0135\times {10}^5}\right)\exp \left(-\frac{0.003475}{T}\right)$

Equations

$$\sum _{\alpha }{\partial }_t\left(\varphi {\rho }_{\alpha }{\chi }_{\alpha }^{\kappa }{S}_{\alpha }\right)+\sum _{\alpha }\nabla \cdot \left({\rho }_{\alpha }{\chi }_{\alpha }^{\kappa }{\mathbf{v}}_{\alpha }\right)=\sum _{\alpha }\nabla \cdot \left(\varphi S_{\alpha }{\mathbf{J}}_{\alpha }^{\kappa }\right)+{\dot{g}}^{\kappa },$$

$${\partial }_t\left(\varphi {\rho }_h{S}_h\right)={\dot{g}}_h.$$

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

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

$${\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 \cdot \left({\rho }_{\alpha }{\mathbf{v}}_{\alpha }{h}_{\alpha }\right)=\nabla \cdot k_{eff}^{th}\nabla T+{\dot{Q}}_h,$$

$$h_{\alpha }={\int }_{T_{ref}}^{T}C{p}_{\alpha }dT,$$

$$u_{\gamma }={\int }_{T_{ref}}^{T}C{v}_{\gamma }dT,$$

$$k_{eff}^{th}=\left(1-\varphi \right)k_s^{th}+\sum _{\beta }\left(\varphi S_{\beta }{k}_{\beta }^{th}\right).$$

$$\sum _{\beta }{S}_{\beta }=1.$$

$$P_g-P_w=P_c(S_w,S_h,\varphi ).$$

$$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.$$

$$\text{Henry’s law,}{z}^{C{H}_4}{\chi }_g^{C{H}_4}{P}_g$$

$$\text{Raoult’s law,}{\chi }_g^{H_{2}O}{P}_g$$

$$\sum _{\kappa }{\chi }_{\alpha }^{\kappa }\le 1\forall \alpha \text{and}\sum _{\kappa }{\chi }_{\alpha }^{\kappa }=1\text{ iff }{S}_{\alpha }>0.$$

$$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 .$$

$${\mathbf{J}}_{\alpha }^{\kappa }=-\tau D_{\alpha }^{\kappa }\left({\rho }_{\alpha }\nabla {\chi }_{\alpha }^{\kappa }\right),$$

$${\dot{g}}^{C{H}_4}=k^r{M}^{CH4}{A}_{rs}\left(P_e-P_g\right)$$

$${\dot{g}}^{H_{2}O}={\dot{g}}^{C{H}_4}{N}_h\frac{M^{H_{2}O}}{M^{C{H}_4}},$$

$${\dot{g}}^{C{H}_4}+{\dot{g}}^{H_{2}O}+{\dot{g}}_h=0.$$

$${\dot{g}}_h<0\text{iff}{P}_g<P_e\text{ and }{S}_h>0,$$

$${\dot{g}}_h>0\text{iff}{P}_g>P_e\text{ and }{S}_g>0,S_w>0.$$

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

$$k^r>0\text{and,}{A}_{rs}$$

$$\text{s.t., }{A}_s$$

$$\text{and, }{\Gamma }_r$$

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

$$P_c=P_{c0}\cdot f_{S_h}^{Pc}\left(S_h\right)$$

$$P_{c0}=p_0{S}_{we}^{-1/\lambda }$$

$$f^{Pc}={\left(1-S_h\right)}^{-\frac{m\lambda -1}{m\lambda }}.$$

$$k_{rw}=S_{we}^{\frac{2+3\lambda }{\lambda }}\text{ and }{k}_{rg}={\left(1-S_{we}\right)}^2\left(1-S_{we}^{\frac{2+\lambda }{\lambda }}\right).$$

$$K=K_0\cdot f^K\left(S_h\right)$$

$$f^K={\left(1-S_h\right)}^{\frac{5m+4}{2m}}$$

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

$$\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)=0,$$

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

Full text

1 Introduction

The motivation for research in methane hydrates is multifarious. Methane hydrates constitute a dominant organic carbon pool in the earth system and an important intermediate ”capacitor” in the global methane budget. Gas hydrates are predominantly formed from biogenic methane that is generated by microbial degradation of organic matter (methanogenesis) in the deep biosphere. This methane is migrating upwards as free gas or methane-rich porewater by advection. This fluid flow is caused by non-steady state sediment compaction (passive margins), compaction of oceanic sediments during subduction (active margins), and dewatering of minerals at elevated temperatures (passive+active margins). Over geological times, the hydrates accumulate close to the bottom simulation reflector (BSR, lower stability limit of gas hydrates) because, the methane flux from below leads to hydrate formation in the gas hydrate stability zone (GHSZ), but 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 anaerobic methane oxidation (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. Our main motivation for modelling the methane hydrate geosystems is to understand this role of gas migration through the GHSZ in the natural carbon cycle. Methane hydrates are also seen as an attractive future energy resource. It is estimated that the total carbon content of methane hydrates is possibly larger than the combined carbon content of all other fossil fuels Pinero2013 ; Burwicz2011 ; Archer2009 . However, there are a number of serious environmental risks associated with the exploitation of gas hydrate reservoirs for the purpose of gas production. Our motivation for modelling the methane hydrate geosystems also extends to the feasibility analysis and risk assessment of various gas production scenarios.

One of the main challenges in modelling the 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/or salinity conditions), 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 disappearence of the free gas phase. The numerical challenges associated with the appearance/disappearance of gaseous or aqueous phases are elaborately discussed in a number of works, like, Class2006 ; Marchand2013 . For methane hydrate models, 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. Secondly, the hydrate and gaseous-aqueous phase transitions manifest at vastly different time scales. For the problems on the geological scales, the rate of methane dissolution is many orders of magnitude higher compared to the rate of hydrate phase change. Typically, the gaseous-aqueous phase transition is modelled as an equilibrium process, while the hydrate phase transition is modelled as a kinetically driven process. Together, both these features of the methane hydrate models compound the numerical challenges of the already complex numerics of phase appearance/disappearance in porous media models.

A number of different numerical methods have been developed to handle the gaseous-aqueous phase transitions in multi-phase multi-component porous media models, e.g., primary variable switching (PVS) schemes WuForsyth2001 ; Class2002 , method of negative saturations Panfilov2014 , method of persistent variables Neumann2013 ; HuangKolditzShao2015 , and non-linear complementary constraints (NCP) approaches Lauser2011 ; Krautle2011 ; GharbiaJaffre2014 ; BuiElman2018 . In the most widely used gas hydrate reservoir simulators, e.g., TOUGH-Hydrate TOUGHHYDRATEv1 , HYRES-C Janicki2011 ; Janicki2014 , STOMP-HYD STOMP , HRS HydrateResSimManual_05 , etc., 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. However, due to the strong coupling and nonlinearities, the phase states in gas hydrate models tend to switch back and forth rapidly, and this often leads to spurious oscillation and a drastic reduction in time step size, in the extreme case, even to a breakdown of the numerical algorithm.

In this article, we present a robust implicit semi-smooth Newton scheme based on an NCP approach for handling phase transitions in methane hydrate models. 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 vapour-liquid-equilibrium (VLE) assumption (e.g.,Helmig1997 ) 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., Facchinei2013 ; Tremolieres2011 ). 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 (See, e.g., HagerWohlmuth2010 and the references therein). We implement our semi-smooth Newton method using an active-set strategy (e.g., Hintermuller2002 ; Hueber2005 ), 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.

In Sec.2, we present our mathematical model and elaborate on the hydrate and the gaseous-aqueous phase transitions. In Sec.3, we introduce our numerical solution scheme for handling these phase transitions. Finally, in Sec.4, we present some numerical examples to validate our numerical model and to show the robustness of our numerical scheme for realistic field scale applications. We also make a comparison of the performance of our semi-smooth Newton scheme with that of a PVS scheme.

2 Mathematical Model

The model is founded on the theory of porous media and considers the reactive transport processes characterizing a typical methane gas hydrate reservoir on the continuum scale. The representative elementary volume (REV) underlying the model is shown in Fig. 1.

The model considers two fluid phases: gaseous and aqueous; and two solid phases: porous granular material (sand or soil) and methane hydrate. The phases are identified with the subscripts $g$, $w$, $s$, and $h$, respectively. We refer to the sand/soil phase as the primary sediment matrix (or simply, the sediment), the hydrate+sediment as the composite sediment matrix, and the fluid and the hydrate phases as the pore-filling phases. The sediment is assumed to be perfectly rigid. The fluid phases are mobile, while the hydrate phase is assumed immobile.

The model considers methane hydrate phase change 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.

A vast majority of methane hydrate geosystems occur in marine settings where the water salinity has strong influences on the thermodynamics and the phase transitions. Therefore, the model also considers the transport of dissolved salts. The model accounts for the miscibility of the fluid phases. Therefore, the model considers that the gaseous phase is comprised of two components: methane and water; while, the aquesous 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, s.t., a single average temperature can be defined over an REV.

In the following, subscript ’$\alpha$’ denotes the fluid phases, while subscript ’$\beta$’ denotes the pore-filling phases, i.e., $\alpha=g,w$, and $\beta=g,w,h$, and the superscript ’$\kappa$’ denotes the components, $\kappa=CH_{4},H_{2}O,c$. The 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$. The fluid phase pressures are denoted with $P_{\alpha}$, the temperature is denoted with $T$, and the porosity is denoted with $\phi$ 111$\phi$ refers to the total porosity, which indicates the void spaces within the primary sediment matrix. This is different from apparent porosity $\phi_{eff}$, which indicates the actual void spaces available for the fluid flow. See Fig.1 for more details..

2.1 Governing equations

2.1.1 Mass, momentum and energy conservation

The transport processes characterizing the gas production from a typical sub-surface methane hydrate reservoir can be 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 HassanizadehGray1979a ; HassanizadehGray1979b ; HassanizadehGray1979c .

Mass balance is considered component-wise for each $\kappa=CH_{4},H_{2}O$,

[TABLE]

where, ${\mathbf{v}_{\alpha}}$ denotes the velocity of the fluid phase $\alpha$ relative to the primary sediment matrix, and ${\bf J}_{\alpha}^{\kappa}$ denotes the diffusive flux of the component $\kappa$ in the phase $\alpha$.

Mass balance for the hydrate phase is given by,

[TABLE]

In Eqn. (1) and (2), the terms $\dot{g}^{CH_{4}}$, $\dot{g}^{H_{2}O}$, and $\dot{g}_{h}$ denote the volumetric source terms resulting from the hydrate phase change.

Mass balance for the dissolved salt is given by,

[TABLE]

where, ${\bf 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 under certain simplifying assumptions (e.g., Helmig1997 ),

[TABLE]

where, $K$ denotes the intrinsic permeability of the composite sediment matrix, i.e., hydrate+sediment matrix, $k_{r\alpha}$ denotes the relative permeabilities, and $\mu_{\alpha}$ the dynamic viscosities.

The primary sediment 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., Helmig1997 ). The energy balance is given by,

[TABLE]

where, $\dot{Q}_{h}$ denotes the heat of hydrate phase change. $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^{th}_{eff}$ is the effective (or lumped) thermal conductivity,

[TABLE]

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,

[TABLE]

The pressures of the fluid phases are related through a capillary pressure $P_{c}$ as,

[TABLE]

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 Sec. 2.2.4.

2.2 Constitutive relations

The $9$ governing equations (1)-(7) consist of the following $25$ unknowns,

[TABLE]

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}^{CH_{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 vapour-liquid equilibrium (VLE), and the Henry’s law and the Raoult’s law are assumed to be valid,

[TABLE]

where, $z^{CH_{4}}$ is the methane gas compressibility, $H^{CH_{4}}_{w}$ is the pressure-corrected Henry’s law solubility constant for methane dissolution in water, and $P^{H_{2}O}_{sat}$ is the saturation vapour 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,

[TABLE]

We can cast the conditions in (10) as a set of Kharush-Kuhn-Tucker complementarity conditions KuhnTucker1951 as,

[TABLE]

2.2.2 Diffusive mass flux

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

[TABLE]

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 \ \sum\limits_{\kappa}{\bf J}_{\alpha}^{\kappa}=0\ hold for all phases $\alpha$. Note that, ${\bf 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 $CH_{4}\cdot N_{h}H_{2}O\rightleftharpoons CH_{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 KimBishnoi1987 , where, the rate of methane and water generated as a result of hydrate phase change are evaluated as,

[TABLE]

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. $M^{\kappa}$ denotes the molar weights, and for methane hydrate, $M_{h}=M^{CH_{4}}+N_{h}M^{H_{2}O}$. Additionally, the following condition holds,

[TABLE]

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

[TABLE]

At $P_{g}=P_{e}$, no reaction can occour, i.e., $\dot{g}_{h}=0$, irrespective of the phase distributions. Also, from Eqns. (13), (14), and (15), it follows that,

[TABLE]

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

[TABLE]

where, $A_{s}$ denotes the specific surface are of the composite sediment matrix, while $A_{0}$ denotes the specific surface are 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 modelled as empirical functions of the form (e.g., Kamath1984 ),

[TABLE]

2.2.4 Hydraulic properties

The capillary pressure $P_{c}$ of the composite sediment matrix is modelled as Gupta2015 ,

[TABLE]

In Eqn.(20), $P_{c0}$ denotes the capillary pressure of the primary sediment 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 BrooksCorey1964 model, where $p_{0}$ is the gas phase entry pressure, $\lambda$ is the soil specific parameter depending on the pore-size distribution, and $S_{we}$ is the normalized aqueous phase saturation, $S_{we}=\dfrac{S_{w}-\left(S_{wr}+S_{gr}\right)}{1-S_{h}-\left(S_{wr}+S_{gr}\right)}$, where, $S_{wr}$ and $S_{gr}$ are the irreducible aqueous and gas phase saturations, respectively.

The relative fluid phase permeabilities are also parameterized following the Brooks-Corey model,

[TABLE]

The intrinsic permeability of the composite sediment matrix is modelled as,

[TABLE]

In Eqn.(22), $K_{0}$ is the intrinsic permeability of the primary sediment 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 SGuptaThesis2016 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\leq 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 Kossel2018 , 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 (Eqn.4) and the constitutive relationships (12)-(19) in the coupled system of Eqns. (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}^{CH_{4}},\chi_{w}^{CH_{4}},\chi_{g}^{H_{2}O},\chi_{w}^{H_{2}O},T\right)\ .$

Eqn. (6) provides an additional relationship for the phase saturations, reducing the number of unknown saturations to $2$. Eqn. (7) provides a relationship for phase pressures, leaving only $1$ pressure unknown. Finally, the Eqn. (8) gives a relationship for $\chi_{\alpha}^{CH_{4}}$ and Eqn. (9) gives a relationship for $\chi_{\alpha}^{H_{2}O}$, thus reducing the unknown mole fractions to $2$. This leaves $7$ primary unknowns which need to be solved for the coupled system which includes $4$ nonlinear second order PDEs (1), (3), and (5), $1$ nonlinear nonhomogeneous first order ODE (2), and $2$ inequality constraints (11).

We choose the following set of primary variables,

[TABLE]

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 gas phase may or may not exist, and along with the hydrate phase saturations, the gas phase saturation and dissolved methane mole fraction are the most important quantities of interest, and therefore, (23) is the most suitable choice.

3 Numerical Solution Strategy

3.1 Space and time discretization of the conservation laws

The Eqns. (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 details of the discretization scheme can be found in SGuptaThesis2016 .

The discretized model can be represented as a system of nonlinear algebraic equations as,

[TABLE]

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 Lauser2011 ,

[TABLE]

which are piecewise linear with respect to the variables $S_{\alpha}$, $\chi_{\alpha}^{\kappa}$. Such functions are commomnly referred as complementary functions or $NCP-$functions in literature. Some examples of other forms of such functions include the minimum function and Fischer-Burmeister function (see Facchinei2013 ; ChenChenKanzow2000 ; Fischer1995a ; Fischer1995b ; Fischer1992 ).

The complementarity constraints (11) and their equivalent form (25) are local in nature, and must hold cell-wise as,

[TABLE]

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

[TABLE]

The active sets $\mathcal{A}_{\alpha}$ corresponds to the cells where phase $\alpha$ is present, while the inactive sets $\mathcal{I}_{\alpha}$ correspond to the calls where phase $\alpha$ is absent.

Using relationships (6), (8), and (9) in Eqn. (26), we get the following system of non-differentiable equations,

[TABLE]

where, ${\Pi_{g}}:=\dfrac{H_{w}^{CH_{4}}}{z^{CH_{4}}P_{g}}$ and ${\Pi_{w}}:=\dfrac{P_{g}}{P_{sat}^{H_{2}O}}$.

3.3 Semismooth Newton scheme

The system of equations (28)-(3.2) 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 HagerWohlmuth2010 . The classical Newton method is valid in all regions where the Eqns. (28)-(3.2) are differentiable, while in other regions where the Eqns. (28)-(3.2) are non-differentiable, the Jacobian can be evaluated by extending the value of the derivatives from the neighbourhood of the non-differentiable regions.

We approximate the Jacobian for our Newton scheme using a central difference method. To approximate the Jacobian for the Eqns. (28)-(3.2), due to their piecewise smooth nature, we use the approximate active/inactive sets $\mathcal{A}_{\alpha}^{(l)}$ and $\mathcal{I}_{\alpha}^{(l)}$ at the $l$-th Newton step to determine the phase wise NCP equations in each cell,

[TABLE]

The approximate active/inactive sets $\mathcal{A}_{\alpha}^{(l)}$ and $\mathcal{I}_{\alpha}^{(l)}$ 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 Sec.3.3 for solving the nonlinear system (24),(28)-(3.2) within the DUNE-PDElab framework Bastian2010 which is based on C++. For solution of the linear system arising from the Newton-linearization, we used a SUPERLU linear solver superLU99_SEQ . For parallel computations, we use 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-in 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 three numerical examples. In Example 1, we validate our numerical model by considering a series of phase transitions involving appearance and disappearance of the gas phase, and comparing the solution with that of a PVS scheme. In Example 2, we simulate the sedimentation driven gas migration through the GHSZ in the highly dynamic geological setting of the Black Sea, and compare the performance of our semi-smooth Newton scheme against a PVS scheme to show the robustness of our numerical scheme. Finally, in Example 3, we simulate a gas production scenario, also based on the geological setting of the Black Sea, where we consider that the hydrate phase is randomly distributed within the hydrate layer with saturations ranging between $0-0.6$, and show that our numerical scheme can robustly handle phase transitions even when the phase distributions and the permeability and porosity profiles are highly complex with large variations.

4.1 Example 1: Model validation

In this example, we verify our numerical scheme by simulating a series of phase transitions over time and comparing the solution with that of a PVS scheme. We start with zero free gas in the domain, and simulate, by manipulating the system pressure, the appearance of the gas phase as a result of hydrate dissociation, followed by the disappearance of the gas phase due to a combination of hydrate reformation and methane dissolution.

4.1.1 Problem setting

We consider a $1m\times 1m$ domain $\Omega$ discretized into $50\times 50$ finite volume cells, and denote the boundary of this domain by $\partial\Omega$. At $t=0$, only the hydrate and the aqueous phases are present in the domain. The gas phase is not present (i.e., $\left.S_{g}\right|_{t=0}=0$), and the aqueous phase contains no dissolved methane (i.e., $\left.\chi_{w}^{CH_{4}}\right|_{t=0}=0$). The hydrate phase is uniformly distributed throughout the domain and has an initial saturation of $\left.S_{h}\right|_{t=0}=0.3$. The initial concentration of the dissolved salt is $\left.\chi_{w}^{c}\right|_{t=0}=5.5$ mmol/mol of water, and the initial temperature in the domain is $\left.T\right|_{t=0}=4^{o}C$. At the initial temperature, pressure, and salinity conditions, the hydrate equilibrium pressure is $\left.P_{eq}\right|_{t=0}=3.4\text{ MPa}$, which is higher than the initial gas pressure $\left.P_{g}\right|_{t=0}=2.0848\text{ MPa}$. So, the hydrate is in an unstable state.

For all $t>0$, the temperature at the boundary is maintained at the initial value, i.e., $\left.T\right|_{\partial\Omega}=4^{o}C$, and a zero-gas-flux condition is prescribed at the boundary, i.e., $\left.\mathbf{v}_{g}\right|_{\partial\Omega}=0$. The phase transitions are triggered by manipulating the boundary conditions for $P_{w}$.

The initial and the boundary conditions are summarized in Table 1. The relevant material properties are listed in Table 4. The problem setting is chosen such, that the spatial variations across the domain are negligible, and the focus of the problem remains on the phase transitions.

4.1.2 Phase transitions

The manipulation of the boundary conditions for $P_{w}$ corresponds to the following four stages (refer Fig. 2(a)):

For the period $0<t<200\text{ hrs}$, $\left.P_{w}\right|_{\partial\Omega}$ is held constant at $2$ MPa. Since hydrate is unstable at this pressure, it will dissociate to produce $CH_{4}$, which will dissolve into the porewater as long as $\chi_{w}^{CH_{4}}<\chi_{w,sat}^{CH_{4}}\left(P_{w},T,X_{c}\right)$, where, $\chi_{w,sat}^{CH_{4}}$ refers to the solubility of methane in the aqueous phase. If the solubility is reached, i.e., $\chi_{w}^{CH_{4}}=\chi_{w,sat}^{CH_{4}}\left(P_{w},T,X_{c}\right)$, the gas phase will appear. 2. 2.

Next, for the period $200\text{ hrs}\leq t<350\text{ hrs}$, a zero-water-flux condition is prescribed, i.e., $\left.\mathbf{v}_{w}\right|_{\partial\Omega}=0$, s.t., the domain is now fully closed. Under these conditions, hydrate will continue to dissociate and pore-pressures will rise until a state of equilibrium is reached s.t., $P_{g}=P_{eq}$. 3. 3.

At $t=350\text{ hrs}$, $\left.P_{w}\right|_{\partial\Omega}$ is instantaneously stepped-up to $5$ MPa, and this pressure is maintained for the period $350\text{ hrs}<t<450\text{ hrs}$. Due to a step increase in the pressure and following the VLE assumption, the gaseous methane will dissolve instantaneously into the porewater, and the gas phase may or may not disappear, depending on how high the new solubility is. If the gas phase does not fully disappear at $t=350\text{ hrs}$, the remaining methane in the gas phase will react with the porewater to form hydrate until the gas phase vanishes. 4. 4.

Finally, for $t\geq 450\text{ hrs}$, $\left.P_{w}\right|_{\partial\Omega}$ is linearly ramped up at a rate of $10$ Pa/s. If the gas phase is still present at $t=450\text{ hrs}$, methane dissolution as well as hydrate formation will continue until the gas phase vanishes and only the aqueous and the hydrate phases remain.

4.1.3 Numerical simulation and results

The numerical simulation was run until $t_{end}=600\text{ hrs}$. The maximum time step size was chosen as $dt_{max}=3600$ s. An adaptive time-stepping strategy was used where the time step size is controlled heuristically based on the number of Newton iterations per time integration step. If the number of Newton iterations is more than $\ell_{h}$, the time step size for the next time integration step is decreased by $25\%$, whereas, if the number of Newton iterations is less than $\ell_{l}$, the time step size is increased by $10\%$. Between $\ell_{l}$ and $\ell_{h}$ Newton iterations, $dt$ is not changed. The choice of $\ell_{l}$ and $\ell_{h}$ depends on the problem setting and, as a rule of thumb, in our numerical scheme we consider $\ell_{l}\geq 4$ and $\ell_{h}=\ell_{l}+4$. In this example, we chose $\ell_{l}=4$ and $\ell_{h}=8$.

The numerical results are shown in Fig. 2, where the $P_{w}$, $S_{g}$, $\chi_{w}^{CH_{4}}$ profiles and the phase state at the point $\left(0.5m,0.5m\right)$ are plotted over time. For the phase state, a value of [math] indicates that the gas phase is present, while a value of $1$ indicates that the gas phase is absent. The results show that in the first stage, the gas phase appears at $t=52\text{ hrs}$. Between $52\text{ hrs}<t<200\text{ hrs}$, $S_{g}$ increases as the hydrate continues to dissociate (Fig. 2(b)). In the second stage, an equilibrium state is achieved at $t\approx 330\text{ hrs}$, s.t., between $330\text{ hrs}<t<350\text{ hrs}$ no further gas dissolution and hydrate phase changes occur (Figs. 2(b),2(c)). In the third stage, at $t=350\text{ hrs}$, solubility of methane is too small for the gas phase to vanish. Between $350\text{ hrs}<t<450\text{ hrs}$, the pressure is constant at $5$ MPa, and a steady state is reached for the dissolved methane, i.e. the rate of gas dissolution equals the rate of hydrate formation, s.t., $S_{g}$ decreases as hydrate formation continues (Fig. 2(b)), while $\chi_{w}^{CH_{4}}$ remains constant (Fig. 2(c)). At $t=450\text{ hrs}$, gaseous methane is still present. As the pressure is ramped up in the fourth stage, both hydrate formation as well as gas dissolution continue until the gas phase finally disappears at $t=482\text{ hrs}$.

In order to ensure that our implementation of the NCP approach for the phase transitions is correct, we compare our results with the more common primary variable switching (PVS) approach of Class2002 . 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 the schemes, we used the same discretization scheme (Sec. 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 compare the solution of our NCP scheme with the PVS scheme, we performed all numerical simulations only in a seuential mode.

We can see in Fig. 2 that both the NCP and the PVS schemes are in agreement over the predicted sequence of the $hydrate\leftrightharpoons gaseous\ methane\leftrightharpoons dissolved\ methane$ phase transitions.

4.2 Example 2: Gas migration through gas-hydrate stability zone (GHSZ)

In this example, we simulate 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 aim 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 compare the performance of our semi-smooth Newton scheme with that of a primary variable switching scheme.

4.2.1 Problem setting

The geological setting for this problem is 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 Winguth2000 . For our problem, we have chosen a buried channel-levee (BCL) complex (blue color, Fig. 3) west of the Viteaz canyon, the main Danube paleo channel, which has buried the BCL over the past $75$ ka (green color, Fig. 3). The BCL is believed to have deposited its levees essentially in two main events correlating to oxygen isotope stages $8$ and $6$ Winguth2000 , i.e. between $320$ ka and $75$ ka BP (brown color, Fig. 3; Zander2017 ). 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 DegensRoss1974 . For the past $500$ ka, this drill core identifies five marine stages that are interrupted by four limnic stages, i.e. 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 Marmara and Mediterranean Sea.

In our problem, we are 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, i.e. in the blue sediments (Fig. 3).

Hence, we base the initial setting $\left(\text{at time }t=t_{0}=0\right)$ on the paleo conditions existing at $300$ ka BP. We choose an arbitrary $1D$ segment $A-B$ located in the eastern levee as our computational domain (See Fig.3). 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 assume 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^{0}$C, corresponding to glacial conditions in the Black Sea Zander2017 . We assume 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^{0}C/km$. The initial conditions for all the primary variables are listed in Table 2. Based on the initial pressure, temperature, and salinity conditions, we can estimate the location of the base of the GHSZ (bGHSZ), i.e., 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}\geq P_{eq}$, and unstable below. (See Fig.4-a.) For this setting, the initial bGHSZ is located $400$m below point A, and we assume that a hydrate layer of $80$m thickess and $30\%$ peak saturation is located just above this initial bGHSZ.

For the sake of simplicity, we assume that the deposition of the brown, and green sediment layers occurs continuously over $300,000$ years at a constant sedimentation rate $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. Fig.4-b 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_{\text{PSF-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 assume that $\Delta t$ is 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 are are fixed at $\left.P_{w}\right|_{z_{\text{PSF-n}},t_{n}}=P_{sf}$ and $\left.T\right|_{z_{\text{PSF-n}},t_{n}}=T_{sf}$, respectively. Note that here we ignore 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, increase over time, which in turn shifts the base of the GHSZ upwards. (Refer to Table 2 for a list of the boundary conditions, and Table 4 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 bel thow 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 stauration 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 single phase model and two-phase model with respect to the $CH_{4}-H_{2}O$ system, especially when the gas phase appears in the domain for the first time.

4.2.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 $dt_{max}=10$ years, and the time step size $dt$ was controlled adaptively using the heuristic time stepping strategy described in Sec. 4.1.3 with $\ell_{l}=8$ and $\ell_{h}=12$.

We performed the numerical simulations with our semi-smooth Newton (NCP) scheme, and for comparison, also with a primary variable switching (PVS) scheme as discussed in Example 1 (Sec. 4.1).

In Fig.7(a), we can see that the NCP scheme took roughly $120$ CPU-hours to solve the problem upto $t=300,000$ years, whereas, due to drastic reduction in time step size, the PVS scheme could reach only upto $t=135,000$ years in twice as many CPU-hours, which is despite the fact that the PVS scheme needs less time per calculation due to fewer degrees of freedom compared to the NCP scheme.

In Fig.5, the snapshots of $S_{h}$, $S_{g}$, $\chi_{w}^{CH_{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 upto which the PVS scheme could solve in $240$ CPU-hours (at which point the simulation was aborted due to large run time). The results of the PVS and the NCP schemes show a very good match, 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 rises 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 Fig.6, this process of hydrate dissociation $\rightarrow$ gas migration $\rightarrow$ hydrate reformation is shown in greater detail by zooming in on the time axis between $60,000$ years $\leq t\leq$ $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 Fig.6(c). This new hydrate layer continues to grow using the free gas supplied by the dissociation of the old hydrate layer, as shown in Figs.6c,d,e.

In Fig.7(b), 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 the schemes. However, the reduction of $dt$ for NCP scheme is not as severe as that for 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.3 Example 3: Gas production through depressurization

In this example, we numerically simulate a gas production scenario where the gas hydrate reservoir is destabilized through depressurization. We consider a single-well configuration, and base the model parameters, material properties, and initial conditions within the reservoir on the geological setting of the Black Sea, similar to Example 2 (Sec.4.2). The objective of this example is 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.3.1 Problem setting

We consider a $2D$ axisymmetric domain, $\Omega$, having dimensions $1000$m $\times$ $1000$m, as shown in Fig. 8. The sea floor, $\partial\Omega_{sf}$, is prescribed at $z=0$m. The depressurizarion well, $\partial\Omega_{well}$, is located at $r=0$m, [math]m $\geq z\geq$ $-400$m. The bottom water temperature at the sea floor is $T_{sf}=4^{0}$C, and the hydrostatic pressure at the sea floor is $P_{sf}=15$ MPa. We assume that the absolute intrinsic permeability and the total porosity of the primary soil skeleton are $K=10^{-13}\text{ m}^{2}$ and $\phi=0.3$, respectively.

At $t=0$, we assume that the domain is fully saturated with saline water, and there is no free gas phase in the domain. Also, the aqueous phase contains no dissolved methane. For the aqueous phase pressure, we consider a hydrostatic pressure gradient, and for the temperature, we prescribe a regional geothermal gradient along the depth, $d_{z}T_{G}=35^{0}$C/km. The bGHSZ is located at $z=-400$m. We consider that an $80$m thick gas hydrate layer, $\Omega_{H}$, exists right above the bGHSZ. To show the robustness of our scheme with respect to complex phase transitions, we prescribe a random distribution of the hydrate phase within this layer, s.t., the hydrate saturation ranges from [math] to $0.6$, and the corresponding absolute permeability ranges from $10^{-13}\text{ m}^{2}$ to $1.6\times 10^{-15}\text{ m}^{2}$.

For $t>0$, a pressure of $\left.P_{w}\right|_{\partial\Omega_{well}}=8$ MPa is prescribed at the production well to simulate gas production through depressurization. At the sea floor, the temperature, pressure, and salinity conditions remain constant and equal to the initial values. At the bottom boundary, $\partial\Omega_{B}$, the regional geothermal gradient is maintained.

The initial and the boundary conditions are listed in Table 3, and the material properties and model parameters are listed in Table 4.

4.3.2 Numerical simulation and results

The computational domain was discretized into a total of $20,276$ quadrilateral elements. The gas production process aws simulated until $t_{end}=360$ days. The maximum time step size was chosen as $dt_{max}=36,000$ sec, and the time step size $dt$ was controlled adaptively using the heuristic strategy discussed in Sec. 4.1.3 with $\ell_{l}=4$ and $\ell_{h}=8$. The simulation was run in parallel on $4$ processing units and required a total of $20$ CPU-hours.

We identify a domain of interest, $\Omega_{I}$, for the gas production process as: [math]m $\leq r\leq$ $250$m, $-150$m $\geq z\geq$ $-550$m. Outside this domain, pressure, temperature, and saturation profiles do not change much. The large size of the domain, however, 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 (QoI) for gas production in gas hydrate reservoirs are $S_{h}$, $S_{g}$, $\chi_{w}^{CH_{4}}$, and the GHSZ. The snapshots of the QoI within $\Omega_{I}$ are plotted in Fig.9 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 is effectively depressurized, and the hydrate layer is fully dissociated within a zone of roughly $15$m around the production well. Due to the relatively high pore pressures, most of the methane is produced from the aqueous phase, and the saturation of the free gas phase remains 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 is the same as described in Sec.4.3.1, except that the hydrate distribution in $\Omega_{H}$ is homogeneous and has a uniform saturation of $0.6$. A snapshop of the QoI in $\Omega_{H}$ for the reference test case is plotted in Fig.10 at time step $t=360$ days. By comparing the behaviour of the numerical solution with the reference test case, we can ensure that the random hydrate distribution has not introduced any artificial numerical artifacts in the numerical solution. In Fig.11, 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 is able to handle the phase transitions quite robustly without significant loss in performance, as indicated by the evolution of the time step sizes.

5 Conclusion

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 modelled as a non-equilibrium based kinetic process, and the phase transitions within the $CH_{4}-H_{2}O$ system were modelled as a VLE process. The inequality conditions resulting from the VLE assumption were cast as 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 Gibson2014 to develop a semi-smooth Newton strategy for solving a Stefan-type problem involving equilibrium based hydrate phase transition.

In Example 1 (Sec.4.1), we validated our semi-smooth Newton scheme against a PVS scheme by simulating a sequence of phase transitions involving appearance and disappearance of the gas phase over time.

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 therfore, can perfom numerical calculations faster. However, in the case of gas hydrate models, this advantage is most often lost because the phase transitions are highly coupled and the PVS scheme requires much smaller time-steps for convergence. We demonstrated this in Example 2 (Sec.4.2) where we simulated gas migration through GHSZ in the highly dynamic geological setting of the Black Sea over paleo time-scales. Notice in Table 4 that, we greatly simplified the problem setting for Example 2 by neglecting the functional dependence of the material properties 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 relatively poorly compared to the semi-smooth Newton scheme. In Example 3 (Sec.4.2), we considered another very important application of methane hydrate models, viz., gas production through depressurization. In this example, we considered a random distribution of the hydrate phase and included strongly nonlinear functional dependencies of the material properties on the local thermodynamic state (see Table 4) in order to show that our semi-smooth Newton scheme can robustly handle even very complex field-scale problems.

Acknowledgements.

This project was funded by the Cluster of Excellence “The Future Ocean”. The Future Ocean is funded within the framework of the Excellence Initiative by the Deutsche Forschungsgemeinschaft (DFG) on behalf of the German federal and state governments.