\documentclass[reqno]{amsart}
\usepackage{graphicx}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
Seventh Mississippi State - UAB Conference on Differential Equations and
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conference 17 (2009),  pp. 1--11.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2009/Conf/17\hfil Methane hydrate]
{Methane hydrate formation and decomposition}

\author[V. Alexiades\hfil EJDE/Conf/17 \hfilneg]
{Vasilios Alexiades} 

\address{Vasilios Alexiades \newline
Department of Mathematics, University of Tennessee,
 Knoxville, TN 37996, USA \newline
 and  Oak Ridge National Laboratory, Oak Ridge TN 37831, USA}
\email{alexiades@utk.edu}

\thanks{Published April 15, 2009.}
\subjclass[2000]{92C45, 35K60, 65M99}
\keywords{Phase change; enthalpy; heat conduction; diffusion; thermochemistry;
 \hfill\break\indent
gas hydrates; finite volume scheme}

\begin{abstract}
Methane hydrates, in arctic permafrost and deep ocean sediments,
store vast amounts of methane, which is the primary constituent of 
natural gas and a potent greenhouse gas.

Methane hydrate is a crystaline solid consisting of methane
surrounded by frozen water molecules. It is stable in a narrow range 
of high pressures and low temperatures. Thus, issues affecting the 
stability of hydrate layers and phase change processes that may 
disturb this stability are of utmost importance.

We present simulations with an isobaric compositional thermal model,
based on conservation laws for species and energy,
which allows composition, temperature, and pressure dependence
of material properties, with thermodynamically consistent treatment
of phase behavior via equations of state.
\end{abstract}

\maketitle
\numberwithin{equation}{section}


\section{Introduction}\label{S:1}

Methane hydrate is ice that burns! (Fig. \ref{fig:1}).
It is a crystalline compound of methane and water, stable under a 
narrow zone of low temperatures and high pressures.
The {\it Structure I hydrate} consists of 46 water molecules per 8 gas 
molecules, forming two small dodecahedral and six large 
tetrakaidecahedral cavities into which methane molecules enter and 
stabilize the structure (Fig. \ref{fig:2}a).  

\begin{figure}[ht]	%%%%%%%%%%%%%%%%%%%%%%% fig 1 a,b
\begin{center}
\includegraphics[width=0.43\textwidth]{fig1a}\quad
\includegraphics[width=0.45\textwidth]{fig1b}
\end{center}
\caption{Ice that burns!} \label{fig:1}
\end{figure}

Its mean radius is about 12 angstroms
and its theoretical composition is 8:46 = 14.8\% CH${}_4$. 
The structure packs a great amount of methane; a liter of {\it pure} 
methane hydrate releases 164 liters of gas!

Methane hydrates are found in sandy sediments at the bottom of oceans 
all around the margins of continents at great depths of 300--3000 m, 
at temperatures of 0--20 ${}^{\circ}$C (Fig. \ref{fig:2}b), 
and under permafrost in arctic regions.
A liter of sediment yields 10 to 50 liters of methane gas.  Thus, 
methane hydrate deposits constitute a potentially enormous natural
gas resource, if the gas could be extracted safely and economically. 
On the other hand, methane being a potent green house gas,
it is feared that global warming may destabilize methane hydrate 
deposits and release vast amounts of methane to the atmosphere.
Thus, issues of formation/dissociation and (thermodynamic) stability 
of hydrates are of utmost concern and need to be studied via 
experiments and modeling.

Primary references on gas hydrates are the books by Makogon 
\cite{makogon97} and by Sloan \cite{sloan98}, while the research
literature is expanding rapidly, see e.g. \cite{wilder08}.
The current state of hydrate research efforts in the USA can be found
in \cite{netl-doe}.

\begin{figure}[ht]	%%% fig 2 a,b
\begin{center}
\includegraphics[width=0.52\textwidth]{fig2a}\quad
\includegraphics[width=0.38\textwidth]{fig2b}
\end{center}
\caption{Methane hydrate crystaline structure (left), and 
temperature-depth stability region (right). From \cite{woodshole}. }
\label{fig:2}
\end{figure}

As a first step in understanding such issues, we present 
a mathematical model of hydrate formation/decomposition under constant
pressure (isobaric) that incorporates compositional and thermal effects,
consistent with the thermochemistry of phases as dictated by the
composition - temperature phase diagram 
of the water--methane binary system $(H_2 O)_{1-x}(CH_4)_x$\,,
with $x$ the mole fraction of methane.

\section{Isobaric Compositional-Thermal Model}\label{S:2}

At a typical pressure of 4.8 MPa, corresponding to a depth of 470 m
(for which the phase diagram is known, \cite{sloan98}),
the water--methane binary $(H_2 O)_{1-x}(CH_4)_x ~$
can form three primary stable phases, (Fig. \ref{fig:phdiag})
which we label as:
\begin{center}
{\bf LG} = {\it Liquid water + methane Gas}, \\
{\bf SL} = {\it Solid hydrate + Liquid water}, \\
{\bf SG} = {\it Solid hydrate + methane Gas}.
\end{center}

At temperature $T_E=6.285^{\circ}$C, the three phases coexist. Above 
$T_E$, water and gas mix at all proportions (except very near $x=0$ 
and $x=1$, shown exaggerated in Fig. \ref{fig:phdiag})).
For $0^{\circ}{\text C} \le T \le T_E$ and $x < x_E=8/54=14.8\%$ 
we find SL (hydrate+water), 
whereas for $x > x_E$ we find SG (hydrate+gas).

\begin{figure}[ht]
  \includegraphics[width=0.7\textwidth]{fig3}
\caption{
Methane hydrate phase diagram (schematic) at 4.8 MPa, 
around the 3-phase coexistence line at $T_E=6.285^{\circ}C$.
}\label{fig:phdiag}
\end{figure}


Employing conservation laws for species and energy,
the $x-T$ phase diagram, and Equations of State (EoS) developed below, 
we seek to determine $x$, $T$, and the mole fractions 
$\lambda^{L}$, $\lambda^{S}$, $\lambda^{G}$ of liquid(L), hydrate(S), 
gas(G), as they evolve in space and time
from given initial and boundary conditions. 

 
\subsection{Conservation laws (for constant $P$, $\rho$)}\label{S:2.1}

Denoting the species as $A=H_2 O$, $B=CH_4$, and their 
molecular weights by $W_A$, $W_B$, the formula molecular weight for
the binary $A_{1-x}B_x$ is $W(x)$=$(1-x) W_A + x W_B$. 
The mass fraction of methane, corresponding to mole fraction $x$, is
$w \equiv w_B$=$x W_B / W(x)$, and that of water is
$w_A$=$(1-x) W_A / W(x)$ since $w_A + w_B$=$1$. 
Conversely, knowing $w$, we can find the mole fraction 
$x \equiv x_B$=$w W_A / [w W_A + (1-w) W_B]$.
Thus, either $x$ or $w$ can be used as composition variable, and 
one can be found from the other. Thermochemistry uses $x$=mole fraction,
but we prefer to express conservation laws in terms of $w$=weight 
fraction, as is commonly done.

Under the assumptions of constant pressure and constant density=$\rho$,
letting $H$=molar enthalpy and $h$=$H / W(x)$=enthalpy (per gram),
conservation of species B (methane) and energy (enthalpy) are expressed 
as:
%
\begin{equation}
 \partial{}_t(\rho w) + \nabla \cdot \vec{J}_B = 0 , \qquad
 \partial{}_t(\rho h) + \nabla \cdot \vec{Q} = 0 ,
\end{equation}\label{eq:1}
with mass flux $\vec{J}_B = - \rho D \nabla w$,
and heat flux $\vec{Q}$ given by
\begin{equation}
  \vec{Q} = - k \nabla T + \overline{h}_A \vec{J}_A + \overline{h}_B \vec{J}_B
 ~= - k \nabla T + \overline{h} \vec{J}_B , \quad
        with \overline{h} = \overline{h}_B - \overline{h}_A.
\end{equation}\label{eq:2}
%
Here $D$=diffusivity of methane in water, 
$k$=thermal conductivity, and 
$\overline{h}_i$=partial enthalpy of species $i=A,B$.
Note that $w_A=1-w$ and $\vec{J}_A = - \vec{J}_B$, so $A$ is also
conserved.

The conservation laws are valid mathematically, in weak sense,
throughout the system, irrespective of phase
(with coefficients appropriate to phase:
$k^{LG}$, $k^{SL}$, $k^{SG}$, $\overline{h}^{LG}$, $\overline{h}^{SL}$, 
$\overline{h}^{SG}$, etc).
%
The conservation laws update $\rho w$ and $\rho h$ to new time,
from which we find $w$ and $h$, and therefore $x$ and $H$ as above.
Then, to find the new phase, temperature, and phase fractions,
we need Equations of State  $H = H(x,T,P)$  consistent with the 
phase diagram (Fig. \ref{fig:phdiag}), which we now develop.


\subsection{Equations of State}\label{S:2.2}

Choose the reference state to be {\it hydrate} at ($x_E$, $T_E (P)$)
with reference enthalpy $H_E$. 
The enthalpy in each phase (LG , SL, SG) can be expressed in terms 
of the pair ($x,T$) by integration, along constant $x$ and constant $T$ 
paths on the phase diagram, of the basic Gibbs relation
%
\begin{equation}
 d H^\alpha = \overline{H}^\alpha dx + C_p^\alpha dT \quad \alpha = \text{LG, SL, SG} ,
\end{equation}\label{eq:3}
%
where 
$\overline{H}^\alpha$=$\overline{H}_B^\alpha - \overline{H}_A^\alpha$, $C_p^\alpha$=heat capacity of phase $\alpha = \text{LG, SL, SG}$.


We obtain the following EoS for $H$ in terms of ($x,T$): \\
%
$\bullet$ In {\bf LG}\quad (for $T \ge T_E$ , $0 \le x \le 1$):
\begin{center} 
$H^{LG}(x,T) = {\bf{\text H^L}} (x) + \int_{T_E}^T {C_p^{LG} (x,\tau)} d\tau $, \\ 
where \quad
${\bf{\text H^L}}(x) := H_E + \Delta H_E^{fus} + \int_{x_E}^x {\overline{H}^{LG} (\xi, T_E)} d \xi $.
\end{center} 
%
$\bullet$ In {\bf SL}\quad (for $T < T_E$, $0 \le x \le x_E$):
\begin{center} 
$H^{SL}(x,T) = {\bf{\text H^S}}(x) - \int_{T}^{T_E} {C_p^{SL} (x,\tau)} d\tau$, \\ 
where \quad
${\bf{\text H^S}}(x) := H_E + \int_{x_E}^x {\overline{H}^{SL} (\xi, T_E)} d \xi$.
\end{center} 
$\bullet$ In {\bf SG}\quad (for $T < T_E$, $x_E \le x \le 1$):
\begin{center} 
$H^{SG}(x,T) = {\bf{\text H^S}}(x) - \int_{T}^{T_E} {C_p^{SG} (x,\tau)} d\tau$, \\ 
where \quad
${\bf{\text H^G}}(x) := H_E + \int_{x_E}^x {\overline{H}^{SG} (\xi, T_E)} d \xi$.
\end{center} 


Since the heat capacity is strictly positive, the dependence of $H$ on
$T$ is monotonic, so $T$ can be found from $H$. Moreover, the quantities
${\bf{\text H^L}}(x), {\bf{\text H^S}}(x), {\bf{\text H^G}}(x)$ 
can be used as switches to distinguish the phases.
Thus, the phases are equally well characterized by the pair ($x,H$)
as follows: \\
%
$\bullet$ If {\bf $H \ge {\bf{\text H^L}}(x)$} then phase is {\bf LG}:

     find $T$ ($\ge T_E$) by solving the EoS 
     $\int_{T_E}^T {C_p^{LG} (x,\tau)} d\tau = H - {\bf{\text H^L}}(x)$ \quad ($\ge 0$) for $T$, 

and then the phase fractions are:
\begin{center}
  $\lambda^{L} = \frac{x^G (T) - x}{x^G (T) - x^L (T)}$, \quad
  $\lambda^{G} = 1-\lambda^L$, \quad
  \quad $\lambda^{S} = 0$.
\end{center}

\noindent
Otherwise, the phase depends on the value of $x$:

\noindent
$\bullet$ For $x < x_E$, if {\bf $H < {\bf{\text H^S}}(x)$}, then phase is {\bf SL} 

     find $T$ ($ < T_E$) by solving the EoS 
     $\int_T^{T_E} {C_p^{SL} (x,\tau)} d\tau = H - {\bf{\text H^S}}(x) \quad (<0)$ for $T$, 

and then the phase fractions are:
\begin{center}
  $\lambda^S = \frac{x - x^L (T)}{x^S (T) - x^L (T)}$, \quad
  $\lambda^L = 1-\lambda^S$, \quad
  $\lambda^G = 0$, \quad if $x < x_E$,
\end{center}

\noindent
$\bullet$ For $x > x_E$, if {\bf $H < {\bf{\text H^G}}(x)$}, then phase is {\bf SG} 

	find $T$ ($ < T_E$) by solving the EoS
     $\int_T^{T_E} {C_p^{SG} (x,\tau)} d\tau = H - {\bf{\text H^G}}(x) \quad (<0)$ for $T$, 

and then the phase fractions are:
\begin{center}
  $\lambda^S = \frac{x^G (T) - x}{x^G (T) - x^S (T)}$, \quad
  $\lambda^G = 1-\lambda^S$   \quad
  $\lambda^L = 0$, \quad if $x > x_E$.
\end{center}

\noindent
$\bullet$ {\bf Otherwise, three-phase L+S+G}, so $T = T_E$,

        and phase fractions are found from the relations:

 $\lambda^L x^L + \lambda^S x^S + \lambda^G x^G = x$, \quad
 $\lambda^L H^L + \lambda^S H^S + \lambda^G H^G = H$, \quad
 $\lambda^L + \lambda^S + \lambda^G = 1$.


The PDEs are discretized by finite volumes, and explicitly in time.
The enthalpy algorithm proceeds as follows:
Knowing the current values of $T$ and phase fractions, 
we evaluate the property values, compute the fluxes $\vec{J}$,
$\vec{Q}$, and update $\rho w$, $\rho h$ from the PDEs. 
From $w$ and $h$ we find $x$, $H$ and the phase indicators
${\bf{\text H^L}}(x)$, ${\bf{\text H^S}}(x)$, ${\bf{\text H^G}}(x)$,
which determine the temperature and phase fractions as just described. 


The model needs the following thermophysical data 
(as functions of $x, T$):
thermal conductivities $k^S$, $k^L$, $k^G$, 
diffusivities (of methane) $D^S$, $D^L$, $D^G$,
partial enthalpies 
$\overline{H}_i^S$, $\overline{H}_i^L$, $\overline{H}_i^G$,
of species i=water, methane, 
heat of fusion (latent heat) $\Delta H^{fus}$, 
and the $x - T$ phase diagram (at the desired pressure).
Some of these quantities were obtained/adapted from the literature
(\cite{makogon97}, \cite{sloan98}, \cite{book93}, \cite{zats-buf}), 
some from phase-equilibria software (\cite{eqtest}), and some
were derived as ideal mixtures of their components.
We do not present them here are they are rather involved and messy.
In addition, simulation of heat transfer in the container walls requires
the density, conductivity, and specific heat of the wall material,
which are constants over our short temperature range of interest
(1 - $25^{\circ}$C).

\section{Numerical Simulations}\label{S:3}


We simulate methane hydrate formation and decomposition in a cylindrical
steel pressure vessel (Fig. \ref{fig:vessel}) 
of inner radius 16 and height 91 cm,
with wall thickness 2.2 cm laterally and 14.6 cm at top and bottom.
The vessel is filled up to height of 64 cm, the rest is water vapor.
In the axially-symmetric simulations, in $(r, z)$ coordinates,
a rather coarse mesh of 32 radial and 128 axial 
nodes was used (roughly 0.5 by 0.5 cm), and 32 by 16 nodes in the void
(plus 4 radial and 2$\times$14 axial nodes in the wall), for a total
of $36 \times 172 = 6192$ control volumes.
This is adequate for the scoping studies presented here as
finer meshes did not make noticeable difference in the results.
A typical simulation takes about 5 minutes 
(on an AMD Opteron 252, 2.6 GHz processor).

\begin{figure}[ht]       %%%%%%%%%%%%%%%%%%%%%%% fig 4 vessel
\begin{center}
 \includegraphics[width=0.3\textwidth]{fig4}
\end{center}
 \caption{Pressure vessel (schematic).} \label{fig:vessel} 
\end{figure}

We present (axi-symmetric) simulations for four scenarios below. In the 
two hydrate formation studies (\S\ref{S:3.1}-\ref{S:3.2}) we start at a 
uniform temperature $T_{init}=25^{\circ}$C and cool the system by 
imposing $T_{bry}=1^{\circ}$C at the outer surface of the vessel.
Conversely, in the two hydrate decomposition studies
(\S\ref{S:3.3}-\ref{S:3.4}) we start at $T_{init}=1^{\circ}$C and heat 
the vessel to $T_{bry}=25^{\circ}$C.
Contour plots for temperature and phase fractions are shown in the
region occupied by the phase-change material (water--methane), 
namely for $0 \le r \le 16$\,cm and $0 \le z \le 64$\,cm 
(with z increasing downwards, representing depth).

 
\subsection{Hydrate Formation}\label{S:3.1}
We start with a water+methane mixture at composition $x = 0.1$ 
(that is, 90\% water - 10\% methane) and
temperature $T_{init}=25^{\circ}$C. According to the phase
diagram the system is in the {\bf LG} phase. We impose 
$T_{bry}=1^{\circ}$C on the outer surface of the vessel. 
Fig. \ref{fig:F-TGLS} displays contours of temperature
and of the gas, liquid, and solid fractions at 6 hours.
Temperature history at three heights (bottom, middle, top, 
at half-radius), 
and radial profiles at 6 and 36 hours are seen in Fig. \ref{fig:Ft-r-T}.
The bottom cools most rapidly and the middle the slowest, as expected.
Complete cooling to $\sim 1^{\circ}$C
occurs in about 84 hours, while the phase fractions relax to essentially
uniform values by 78 hours. The result is 70\% solid (hydrate) and 
30\% liquid (water). 

\begin{figure}[ht]       %%%%%%% fig 5 a,b,c,d  F-TGLS
\begin{center}
\includegraphics[width=0.25\textwidth]{fig5a}\;
\includegraphics[width=0.22\textwidth]{fig5b}\;
\includegraphics[width=0.23\textwidth]{fig5c}\;
\includegraphics[width=0.22\textwidth]{fig5d}
\end{center}
\caption{Formation (see \S\ref{S:3.1}):
Contour plots of $T$ and phase fractions  
$\lambda^G,\lambda^L, \lambda^S$ at 6 hours.} \label{fig:F-TGLS}
\end{figure}


\begin{figure}[ht]       %%%%%% fig 6 a,b  Ft,r-T
\begin{center}
\includegraphics[width=0.45\textwidth]{fig6a}\quad
\includegraphics[width=0.45\textwidth]{fig6b}
\end{center}
\caption{Formation (see \S\ref{S:3.1}):
Temperature history at three heights (left). Radial profiles 
at 6 and 36 hours (right).}
\label{fig:Ft-r-T}
\end{figure}


\subsection{Hydrate Formation with gas bubble}\label{S:3.2}
Now we start again in the {\bf LG} phase
at $x = 0.1$, $T_{init}=25^{\circ}$C and impose $T_{bry}=1^{\circ}$C,
but we assume there is a pre-existing methane gas bubble 
($\lambda^G=1$, $\lambda^L=\lambda^S=0$)  at the center. 
Contour plots of temperature and of the phase fractions at 6 hours
are shown in Fig. \ref{fig:Fbbl-TGLS}. 

\begin{figure}[ht]      %%%%%%% fig 7 a,b,c,d  Fbbl-TGLS
\begin{center}
\includegraphics[width=0.24\textwidth]{fig7a}\;
\includegraphics[width=0.22\textwidth]{fig7b}\;
\includegraphics[width=0.225\textwidth]{fig7c}\;
\includegraphics[width=0.22\textwidth]{fig7d}
\end{center}
\caption{Formation with gas bubble at the center initially (see \S\ref{S:3.2}).
Contour plots of $T$ and $\lambda^G,\lambda^L, \lambda^S$ at 6 hours. 
} \label{fig:Fbbl-TGLS} 
\end{figure}

 
Complete cooling to nearly $1^{\circ}$C again takes about 84 hours,
but now the phase fractions do not relax to their final 
uniform values till much much later.
Even after 240 hours there are non-uniformities, shown 
in Fig. \ref{fig:Fbbl-240h}, revealing higher hydrate fraction around
the pre-existing gas bubble. It persists for a very long time after
temperature gradients have dissipated, due to the very low diffusivity 
of methane.
The final result is again 70\% hydrate and 30\% water.

\begin{figure}[ht]     %%%%%%% fig 8 a,b,c  Fbbl-240h
\begin{center}
\includegraphics[width=0.26\textwidth]{fig8a}\quad
\includegraphics[width=0.27\textwidth]{fig8b}\quad
\includegraphics[width=0.26\textwidth]{fig8c}
\end{center}
\caption{Formation with gas bubble at the center initially (see \S\ref{S:3.2}).
The persisting phase fractions $\lambda^G,\lambda^L, \lambda^S$ at 240 hours. 
} \label{fig:Fbbl-240h} 
\end{figure}

\subsection{Hydrate Decomposition}\label{S:3.3}
Here we start in the {\bf SL} phase with 
hydrate+water at $x = 0.1$, at a cold temperature $T_{init}=1^{\circ}$C 
and impose $T_{bry}=25^{\circ}$C on the outer surface of the vessel. 
Contours of temperature and of the gas, liquid, and solid fractions 
at 24 hours are shown in Fig. \ref{fig:D-TGLS}. 
Complete warming to nearly $25^{\circ}$C takes about 168 hours, 
much longer than the corresponding formation case of \S\ref{S:3.1}. 
The result is 89\% liquid (water) and 11\% gas (methane).

\begin{figure}[ht]      %%%%%%%%% fig 9 a,b,c,d  D-TGLS
\begin{center}
\includegraphics[width=0.24\textwidth]{fig9a}\;
\includegraphics[width=0.21\textwidth]{fig9b}\;
\includegraphics[width=0.22\textwidth]{fig9c}\;
\includegraphics[width=0.21\textwidth]{fig9d}
\end{center}
\caption{Decomposition (\S\ref{S:3.3}): Contour plots of temperature and 
$\lambda^G,\lambda^L, \lambda^S$ at 24 hours. 
} \label{fig:D-TGLS} 
\end{figure}

 
\subsection{Hydrate Decomposition with gas bubble}\label{S:3.4}
In this simulation, we start again in the {\bf SL} phase
at $x = 0.1$, $T_{init}=1^{\circ}$C and impose $T_{bry}=25^{\circ}$C,
but now we assume there is an all-gas bubble 
($\lambda^G=1$, $\lambda^S=\lambda^L=0$)  at the center.
Complete heating to nearly $25^{\circ}$C takes about 160 hours,
but the evolution of the phase fractions is quite
different, shown in Fig. \ref{fig:Dbbl-TGLS} at 24 hours. 


\begin{figure}[ht]      %%%%%%% fig 10 a,b,c,d  Dbbl-TGLS
\begin{center}
\includegraphics[width=0.24\textwidth]{fig10a}\;
\includegraphics[width=0.21\textwidth]{fig10b}\;
\includegraphics[width=0.22\textwidth]{fig10c}\;
\includegraphics[width=0.215\textwidth]{fig10d}
\end{center}
\caption{Decomposition with gas bubble at the center (\S\ref{S:3.4}).
Contour plots of temperature and $\lambda^G,\lambda^L, \lambda^S$ at 24 hours. 
} \label{fig:Dbbl-TGLS} 
\end{figure}


The solid phase disappears by 132 hours, but the liquid and gas 
fractions take a very long time to relax to uniform values, 
seen at 240 hours in Fig. \ref{fig:Dbbl-240h}.
The ultimate result is again 89\% water and 11\% gas.

\begin{figure}[ht]       %%%%%%%%%%%%%%%%%%%%%%% fig 11 a,b  Dbbl-240h
\begin{center}
\includegraphics[width=0.24\textwidth]{fig11a}\quad
\includegraphics[width=0.25\textwidth]{fig11b}
\end{center}
\caption{Decomposition with gas bubble at the center (\S\ref{S:3.4}).
Persisting $\lambda^G$ and $\lambda^L$ at 240 hours. 
} \label{fig:Dbbl-240h} 
\end{figure}


\section{Discussion}\label{S:4}
An isobaric conduction-diffusion model for hydrate formation and 
decomposition was presented, based on species and energy conservation 
laws, and Equations of State consistent with thermochemistry and the 
phase diagram of the water--methane binary system 
$(H_2 O)_{1-x}(CH_4)_x ~$. Given the initial composition $x$, 
temperature $T_{init}$, and imposed temperature $T_{bry}$ on the outer
surface of the vessel, 
the model tracks $x$, $T$, and the phase fractions 
$\lambda^G,\lambda^L, \lambda^S$ of gas (methane), liquid (water), 
solid (hydrate), as they evolve in space and time.

Axially-symmetric numerical simulations in a cylindrical pressure vessel
were carried out to find out basic features of hydrate formation and 
decomposition processes. These include how long the process takes 
(for gradients to dissipate), and the final phase fractions produced.

Sample simulations were presented in \S\ref{S:3} for composition
$x=0.1$. The behavior is similar for other compositions we have tested
($x= 0.15, 0.2, 0.4, 0.5, 0.8, 0.9$) under conditions identical to those
in \S\ref{S:3.1} (for formation) and \S\ref{S:3.3} (for decomposition). 
Pertinent resuls are summarized in the tables below.

Table \ref{Tab:duration} shows how long it takes for a formation or
decomposition process to complete for various compositions.
Decomposition is much slower than formation 
(as already noted in \S\ref{S:3.3}), except in the $x = 0.8$ case 
in which the gas content is very high.

\begin{table}[ht]
\caption{\label{Tab:duration}
Number of hours for each quantity to roughly equilibrate to its final 
value, for formation (\S\ref{S:3.1}) and decomposition (\S\ref{S:3.3})
for various initial compositions $x$.}
\renewcommand{\arraystretch}{1.25}
\begin{tabular}{|c||c|c||c|c|}
\hline
& \multicolumn{2}{c ||}{Formation} & \multicolumn{2}{c|}{Decomposition} \\
\hline
$x$ & $T ^\circ$C & $\lambda^G$,$\lambda^L$,$\lambda^S$ \~(\%)     & $T ^\circ$C & $\lambda^G$,$\lambda^L$,$\lambda^S$ \~(\%) \\
\hline
0.10 	& 84 & 78     & 168 & 132    \\
0.15 	& 66 & 66     & 162 & 132    \\
0.20 	& 66 & 66     & 162 & 132    \\
0.50 	& 96 & 90     & 162 & 132    \\
0.80 	& 170 & 166   & 148 & 132    \\
\hline
\end{tabular}
\end{table}

Finally, in Table \ref{Tab:fractions} we list the initial and final
percentage of phases, for various compositions $x$,
in the case of formation (cooling from $25^{\circ}$C down to 
$1^{\circ}$C), 
and of decomposition (heating from $1^{\circ}$C up to $25^{\circ}$C). 

\begin{table}[ht]
\caption{Initial and final phase fractions (expressed in \%)
}\label{Tab:fractions}
\renewcommand{\arraystretch}{1.25}
\begin{tabular}{|c||c|c|c||c|c|c|}
\hline
 & \multicolumn{3}{c ||}{initial values (\%)} & \multicolumn{3}{c|}{final values (\%)} \\
\cline{2-4} \cline{5-7} 
$x$ & $\lambda^G$ & $\lambda^L$ & $\lambda^S$ & $\lambda^G$ & $\lambda^L$ & $\lambda^S$\\
\hline
\multicolumn{7}{| l |}{Formation} \\
0.10 & 11 & 89 & -   &  0 & 30 & 70 \\
0.15 & 16 & 84 & -   &  1 &  0 & 99 \\
0.20 & 21 & 79 & -   &  7 &  0 & 93 \\
0.50 & 54 & 46 & -   & 42 &  0 & 58 \\
0.80 & 86 & 14 & -   & 77 &  0 & 23 \\
\hline
\multicolumn{7}{| l |}{Decomposition} \\
0.10 &  - & 30 & 70  & 11 & 89 &  0 \\
0.15 &  1 &  - & 99  & 16 & 84 &  0 \\
0.20 &  7 &  - & 93  & 21 & 79 &  0 \\
0.50 & 42 &  - & 58  & 54 & 46 &  0 \\
0.80 & 77 &  - & 23  & 86 & 14 &  0 \\
\hline
\end{tabular}
\end{table}

We observe that for each composition $x$ one process undoes the effect 
of the other, recycling the system back to its original state
(under the ideal conditions assumed here).
E.g. formation for $x=0.2$ starts with $\lambda^G = 21\%$, 
$\lambda^L=79\%$ and produces $\lambda^G=7\%$, $\lambda^S=93\%$.
Conversely, decomposition for $x=0.2$ starts with $\lambda^G=7\%$,
$\lambda^S=93\%$ and produces $\lambda^G=21\%$, $\lambda^L=79\%$.

The maximum amount of hydrate is obtained for composition $x=0.15$,
as dictated by the phase diagram,
because $x=0.15$ is very close to the theoretical ratio 
$x_E = 8:46 \approx 0.148$ of hydrate formation.

The model can be used to study sensitivity of parameters and 
uncertainty quantification, and towards determining unavailable
thermophysical parameters.
For increasing realism, the model can be enhanced by incorporating
effects of salinity, of crystallization and decomposition kinetics, 
of sediments, and ultimately also of microbiological activity as a 
source of methane in sediments.


\begin{thebibliography}{00}	

\bibitem{book93}
V Alexiades and A D Solomon
{\bf Mathematical Modeling of Melting and Freezing Processes},
Hemisphere Publ.Co., 1993.

\bibitem{makogon97}
Yuri F Makogon,
{\bf Hydrates of Hydrocarbons}, Pennwell Books, 1997.

\bibitem{netl-doe}
National Methane Hydrate R\&D Program,
http://www.netl.doe.gov/technologies/oil-gas/futuresupply/methanehydrates/maincontent.htm \ \ (accessed September 2008).

\bibitem{sloan98}
E ~Dendy Sloan, Jr.,
{\bf Clathrate Hydrates of Natural Gases}, 2nd Ed., Marcel Dekker, 1998.

\bibitem{eqtest}
W Wagner, EQTEST (fortran code) for thermodynamic properties of H2O and of CH4, Ruhr-Universität Bochum, 1995.

\bibitem{wilder08}
Joseph W Wilder, et.al.,
"An International Effort to Compare Gas Hydrate Reservoir Simulators"
Proceedings of 6th International Conference on Gas Hydrates (ICGH 2008), Vancouver, CANADA, July 6-10, 2008. \ \ 
http://www.netl.doe.gov/technologies/oil-gas/FutureSupply/MethaneHydrates/MH\_CodeCompare/MH\_CodeCompare.html \ (accessed September 2008).

\bibitem{woodshole}
Woods Hole Science Center http://woodshole.er.usgs.gov/project-pages/hydrates/what.html \ (accessed September 2008).

\bibitem{zats-buf}
O.Y. Zatsepina and B.A. Buffett,
"Phase equilibrium of gas hydrate: implications for the formation of
        hydrate in the deep sea floor",
{\it Geophys. Res. Lett.} {\bf 24}: 1567-1570, 1997

\end{thebibliography}

\end{document}