Why Should We Use Residual Thermodynamics for Calculation of Hydrate Phase Transitions?

: The formation of natural gas hydrates during processing and transport of natural has historically been one of the motivations for research on hydrates. In recent years, there has been much focus on the use of hydrate as a phase for compact transport of natural gas, as well as many other applications such as desalination of seawater and the use of hydrate phase in heat pumps. The huge amounts of energy in the form of hydrates distributed in various ways in sediments is a hot topic many places around the world. Common to all these situations of hydrates in nature or industry is that temperature and pressure are both defined. Mathematically, this does not balance the number of independent variables minus conservation of mass and minus equilibrium conditions. There is a need for thermodynamic models for hydrates that can be used for non ‐ equilibrium systems and hydrate formation from different phase, as well as different routes for hydrate dissociation. In this work we first discuss a residual thermodynamic model scheme with the more commonly used reference method for pressure temperature stability limits. However, the residual thermodynamic method stretches far beyond that to other routes for hydrate formation, such as hydrate formation from dissolved hydrate formers. More important, the residual thermodynamic method can be utilized for many thermodynamic properties involved in real hydrate systems. Consistent free energies and enthalpies are only two of these properties. In non ‐ equilibrium systems, a consistent thermodynamic reference system (ideal gas) makes it easier to evaluate most likely distribution of phases and compositions. Δ T below hydrate stability T for a fixed pressure. The


Introduction
The problems of hydrate formation in pipelines during transport of hydrocarbons and other hydrate forming components is as old as the modern oil industry itself. The need for calculations of hydrate formation conditions in order to design appropriate measures to counteract problems of pipeline blockings is a continuous effort. During the last three decades there has been a substantial increase in the interest of natural gas hydrates as an energy source, which requires calculation of phase transition conditions and phase transition kinetics. However, also other sides of natural gas hydrates motivate the developments of better and more complete tools for calculations of hydrate phase transitions. Hydrate exposed to inflow of seawater through fracture systems leads to leakage fluxes of methane to the oceans and potentially also to air. All these dynamics processes may also lead to geo mechanical instability and landslides.
Calculation of hydrate phase transitions has a long history. Many early strategies were based on analogies to calculations of gas/Liquid distributions in hydrocarbon systems. Gas-liquid mole-fraction ratios and K-values are used along with a mass balance to calculate distributions of oil and gas, as well as composition of the phases. In the early days before computers became commercially available, K-value charts for various components were developed. Similar K-value charts were also developed for various hydrate formers in a similar analogy for hydrate equilibrium calculations. It is far beyond the scope of this work to discuss the very old strategies for hydrate equilibrium calculations and for this we refer to some history in Koh & Sloan [1]. Using a semi-grand canonical ensemble, van der Waal & Platteeuw [2] derived a Langmuir-type adsorption theory for hydrate which generated various ways to treat the hydrate phase in a more modern fashion, using equations of state to describe the impact of hydrate formers on thermodynamic equilibrium. One version of the modern way to model hydrate equilibrium was based on the use of a reference hydrate, and mostly credited to Parrish and Prausnitz [1,3]. For review of other methods and more details the reader is directed to other literature such as [1].
In this work we only focus on two different approaches. The first method is what we can call a reference approach since it utilizes a reference state and differences between a pure water phase and empty hydrate of either structure I, II or H. A second method uses residual thermodynamics for all components in all phases, including hydrate.

Motivation and Overview
Frequently independent thermodynamic variables are often used to evaluate energy processes. One typical example is evaluation of risk for hydrate formation, which is frequently discussed in terms of pressure and temperature stability limits. This projection of the whole stability regime, which also includes concentrations in all co-existing phases do not tell directly if the free energy change needed to create hydrate. Moreover, it does tell anything about how the released heat of hydrate formation and dissociation, is a multi-phase problem in which Gibbs distributed away from the formed hydrate. This is just one example that tells us that we need a system for analyzing hydrates which is based on thermodynamic responses rather than independent thermodynamic variables. Practically this means that we need to develop model systems which analyze hydrate formation and hydrate dynamics based on free energy changes for the variety of possible (combined first and second law) as responses to changes in temperature, pressure and concentrations. Any phase transition determine phase distributions under constraints of dynamics is implicitly coupled to related mass and heat transport. Pressure and temperature dynamics and heat transport dynamics. We therefore also need a consistent rote to calculations of enthalpies, as the first law response to independent variables like temperature, pressures and concentrations. Similar examples related to hydrates in porous media are just a two-dimensional projection of all independent thermodynamic variables. Concentrations are also discussed in this study.
Another important motivation for this work is that hydrates in porous media can never reach true thermodynamic equilibrium. As we demonstrate in this study there is a lower limit of all hydrate formers and former in surrounding water in all co-existing phases are additional independent thermodynamic variables. One which is needed in order to keep hydrate stable. There is even a lower limit of the objectives water in gas needed in order to prevent the hydrate from sublimation. Overall, there are far too many active phases of this work is to provide a thermodynamic overview of the relevance to hydrate formation and hydrate stability to fulfill the balance between number of some independent variables on one side and conservation laws plus condition of the phase transitions which are often omitted. thermodynamic equilibrium on the other side. When temperature and pressure are both defined in a reservoir or in a pipeline two independent variables are already fixed. Even with only one hydrate former plus water there are three phases when hydrate forms. That leaves only one independent thermodynamic variable and two are defined. In reality, as also discussed here, hydrate systems in porous media is more mathematically over determined compare tow two fixed independent thermodynamic variables. This is just one example that tells us that we need a new thermodynamic toolbox which is able to calculate all the hydrate phase transitions of significance for hydrates in sediments or hydrates forming during transport of hydrate formers in a pipeline. This also involves the need for calculation of enthalpies of hydrate phase transitions. Moreover, since hydrate formed from different phases has different stability there is need for a toolbox which also calculated free energies of the various co-existing phases. This is a second objective of this work.
Offshore methane hydrate reservoirs are always in a dynamic state. This implies that fracture systems from below bring in hydrocarbons that lead to formation of hydrate with groundwater. At the same time seawater is leaking into the hydrate filled sediments through fracture systems. As will be discussed in more detail in this study this leads dissociation of hydrates because the seawater content of methane is almost zero. A third objective of this work is to shed more light on important hydrate stability limits, beyond the temperature pressure projection of the stability limits There are many thermodynamic packages for calculating pressure and temperature stability limits. These are based on very old calculation routes from around 1970. There are many drawbacks related to these old approaches. And a fundamental limitation is that thermodynamic properties like chemical potentials are empirically fitted. Practically these packages only calculate hydrate formation from a separate hydrate former phase and liquid water or ice. It is not What is new here is a complete concept for calculating hydrate stability limits in various projections. Not only in temperature and pressure, but also for hydrate formation from dissolved hydrate formers and hydrate dissociation towards water under saturated with hydrate former. These types of hydrate phase transitions are critical in the hydrate dynamics related to fracture systems that connects offshore hydrate bearing formations in contact with seawater. However, the most important is a goal of this study to do an extensive comparison of the residual scheme and the reference model. However, it is actually fairly simple to rewrite programs based on the reference over to a residual complete and consistent thermodynamic model system. A fourth objective is therefore to illustrate the thermodynamic similarities, and hopefully, to illustrate what changes are needed to reformulate reference schemes over to residual thermodynamic schemes. That can easily be implemented into hydrate reservoir simulators or hydrate risk evaluation software. This can make substantial steps forward in evaluation of hydrate production scenarios.
The next section gives an overview of the residual thermodynamic concept, along with a brief description of what we call the reference method. The main purpose of this section is to point of advantages and drawbacks of the two different schemes, and also provide a platform for residual thermodynamic analysis of other routes to hydrate dissociation and formation.
In Section 3, we show some examples for hydrate stability limits as based on the residual thermodynamic scheme and on the reference method. Since we do not have any code for the reference method we have used software which is publicly available and documented in many other publications from other groups. Another new element in this section is the estimation of a hydrate curve for CO2 which also includes the impact of a CO2 phase transition which is frequently overlooked and smoothened out. Section 4 focuses on hydrate formation from dissolved hydrate formers as well as the dissociation limits for hydrate in presence of water which in under saturated with hydrate formers. These phase transitions are very important in analysis of hydrate dynamics in sediments, but also during hydrate formation and dissociation in a multiphase flow line containing hydrocarbon oil and gas phases and separate liquid water phase. Hydrate can nucleate and form towards rusty pipeline walls and on water-gas interface (and potentially water-liquid interface). During turbulent flow, these hydrates can dissociate again when exposed to water lean on hydrate formers.
Hydrate dissociation needs two conditions to be fulfilled. The free energy change has to be large enough to efficiently release water and hydrate former from the hydrate. During formation of hydrate there is a barrier related to the push work needed for make space for the new phase. During hydrate dissociation the guest molecules have to cross a thin interface of structured water and correspondingly low diffusivity. In addition to this "penalty" of slow mass transport and need for a significant free energy difference the necessary heat must be supplied. Reducing pressure to below temperature and pressure stability is one method for producing hydrate. The questions are; are the free energy differences sufficient and how is heat supplied. Is there sufficient heat supply? Estimation of consistent enthalpies is crucial in hydrate production. Moreover, the calculations need to be consistent with free energy calculation for the phase transition changes in order to give the correct entropy generation. In Section 5, we discuss thermodynamic models for Gibbs free energy and enthalpy derived from the residual thermodynamic concept.
The study is completed with a discussion in Section 6, followed by our conclusions.

Thermodynamic Models in Residual Thermodynamics Model and the Reference Models Method
In a thermodynamic description we will use m k  as symbol for chemical potential for component k in a phase m. Within the limitations of this work m will be water (ice or liquid), hydrate, gas (hydrate former phase as gas, liquid or supercritical) and adsorbed. Index k will be H2O and any other component that distributes over the phases m. This also include possible thermodynamic inhibitors. Fugacity m k f is defined for each component as: R is the universal gas constant, T is temperature and P is pressure. Equation (1) Equations (1) In the original derivation by van der Waal and Platteeuw [2] the water lattice were assumed to be rigid while a later derivation (Kvamme & Tanaka [4]) permitted movement of the water molecules in the lattice. This latter approach made it possible to investigate the effect of guest movements on the water lattice by using a different evaluation for the impact of the guest movements. The treatment of guest molecule movements in the cavity as a harmonic oscillating spring, from minimum energy state in a molecular dynamics study [4] provides insight into hydrate destabilization effects due to size and mass. Chemical potential for molecule type i in cavity type k. We will assume that small and large cavities are at equilibrium so that: For a system at equilibrium then the chemical potential for a guest molecule in a cavity is equal to the chemical potential for the same molecule in the equilibrium phase. Δgk1 is the Gibbs free energy change for inclusion for guest molecule I in a cavity of type k.
The most classical example is a hydrate former phase (gas, liquid, supercritical) in contact with liquid (or ice) water that form a hydrate. For these three phases there are 12 independent thermodynamic variable, 3 conservation laws and 8 conditions of equilibrium. As is trivially known we can then fix one independent thermodynamic variable, commonly T or P. For this particular case equation of ququilibrium we have: For the equilibrium between the gas and the hydrate.
In the classical formulation of van der Waal & Platteeuw [3], an alternative formulation for (9) for a rigid water lattice is: The Langmuir constant ( ) ki C T for a molecule i in cavity k and given below as Equation (12).
For a molecule like methane the results from (9) and (11) are almost the same while smaller molecules such as N2 are better represented by (11). For a large molecules likw CO2 the difference in impact on water hydrate chemical potential, Equation (6), is one kJ/mole since the movements of CO2 interferes with some water lattice librations. In the simplest case of a monoatomic spherical guest molecules the Langmuir constant is a simple integral over the Boltzmann factors of interaction energies between the guest molecule and surrounding waters. (12) For nonlinear multi-atomic representations of guest molecules the integration will involve rorational degrees of freedom. Guest-guest interactions between guest in different are also significant [5]. Polar guest molecules such as H2S will also get extra stabilization from coulumbic interactions between the partial charges in H2S and water molecules in the cavity lattice [6]. For CO2, on the other hand, the quadropole moment will result in a destabilization effect [6]. However, for now Equation (12) serves as sufficient illustration. Various simplifications of (12) can be found in [2].
The most common guest/water interaction model in present versions hydrate equilibrium codes based of the reference method is based on a spherically smeared out version of the Kihara potential for interactions between a water and a guest. The Kihara potential can be expressed as: 12 6 ( ) i and j are molecular indexes while ij ij r a  is the closest distance between the two molecules. ij  is a molecular diameter and ij  is a well-depth. For aij equal to zero (13) reduces to the Lennard-Jones 12-6 potential which we and many others have utilized in various studies. See for instance references [4][5][6]. A summation of pariwise interactions in Equation (12) is possible and integration can be conducted efficiently using a Monte Carlo approach [5,6], but it is more common to use an integrated smeared interaction version in which the average water/guest interaction are smeared out over the surface of a speheriaclly smoothed cavity radius R with z being the number of waters represented in this spehical shell. Z is therefore 20 for small cavity and 24 for large cavity. The details of this integration to reach at the spherically smoothed potential is far too extensive to include here. See reference [2] for more details and further references. The final results is for each specific cavity k is: The sperically symmetric integration version of (12) can then be expressed as: Some Kihara parameters for the smoother cavity approach are listed below in Table 1. These are of course fitted also with specific fitted parameters when the reference approach is used. As such the Kihara papameters in Table 1 should be used with reference parameters published from the same research groups. List of various published reference properties are listed in Table 1 below. Cavity raidii published and coordination numbers published by various research groups are listed in Table  2 below. In order to complete the rquilibrium calculation for hydrate formation between gas, liquid water and hydrate the symmetric excess formulation of water chemical potential is:  is the activity coefficient of the liquid water as function dissolved hydrate formers as well as additives like methanol and salt. One approach for solving the equilibrium for water is based on residual thermodynamics also for hydrate. For a well defined activity of water accirding to impacts of solutes the solution of Equation (17) is feasible because liquid water chemical potential as well as empty hydrate chemical potential is known from molecular dynamics simulations and verified in many publication. Some recent examples are [22][23][24][25][26][27][28][29][30].
With known gas composition and a model for the gas fugacity coefficient, we have utilized the SRK [31] equation of state in anumber of recent publications. Equations (6) and (17) in (18) can be solved for T if a pressure is given or alternatively for P when temperature is given.
In the absence of data for liquid water chemical potential and water chemical potential for water in empty clathrate of either structure of hydrate Equations (6) and (17) in (18) can be reformulated to: In which either (9) or (11) can be utilized to calculate the cavity partition functions for small and large cavity fillings of the varous guest molecules in the system. Equation (19) is hereafter denited as the reference method. Direct iterative solution of equation (18) using chemical potentials for for pure liquid water (or ice) as well as empty clatharet water chemical potential in (6) from Kvamme & Tanaka [4] is now the residual thermodynamic method. The chemical potential difference in (19) is typically fitted towatds experimental data through the following parameters: H O H T P  is the enthalpy difference between liquid water enthalpy and empty hydrate water enthalpy.
is the specific heat capacity difference between liquid water and empty hydrate for the sppecific structure in consideration.
Liquid water density does not change much over the limited range of liquid water temperatures for hydrate stability. There is a slight temperature dependency in the hydrate lattice constant [32], but not substantial so a constant (22) is fair enough, as also indicated in the equation.
In summary the reference approach needs fitted values for (20), Altogether 5 parameters that needs to be fitted.
Some selected values from open literature for the parameters discussed above is listed in Table  3 below. There may be many more since the various groups using this method may not always publish their fitted values.   [3,13,33,34] In addition to the fitting of fundamental thermodynamic properties the interaction energies between water and guest molecules typically involves fitting of three parameters in a Kihara type of potential for each guest molecule in each type of cavity for the integral in Equation (12). These integrations are normally conducted over a spherically smoothened cavity. See for instance Sloan's book [2]. While the small cavity of structure I is symmetric the large cavirt in structure I is asymmetric and on average non-spherical due to the two hexagonal faces.

Hydrate Stability Limits in the Pressure-Temperature Projection of Independent Thermodynamic Variable
For one hydrate former and liquid water distributed over 3 phases the number of independent thermodynamic variables are 12 and the sum of conservation laws and conditions of equilibrium is 11. Equilibrium is therefore only possible if one thermodynamic variable is defined. For given temperatures we can therefore solve conditions of equilibrium according to (18) using either the residual thermodynamic scheme or the reference scheme. For the latter alternative we could only find CSMHYD [39] as an open source to compare with, along with experimental data. A comparison is plotted in Figure 1 for CH4 hydrate. A comparison for CO2 hydrate is plotted in Figure 2. A comparison between calculated stability limits for a mixture of CO2 and CH4 is plotted in figure 3 and compared to experimental data from open literature. Note that CSMHYD do not estimate the phase transition over to more dense CO2 phase.   Hydrate formed from CO2 gas and liquid water. Calculated quilibrium curve for 3.46% CH4 + 96% CO2 hydrate from this work (solid) as compared with CSMHYD (dashed line) [39] and experimental data from Fan and Guo (1999) [45] (circles).
Even if another hydrate former is added so that Gibbs phase rule is achieved it does not mean that the system can reach equilibrium. The reason is three-fold: (1) More than one hydrate phase forms due to formation from separate hydrate former and water pluss hydrate forming from dissolved hydrate former in water or hydrate former adsorbed on minerals; (2) Even for hydrate forming from a separate hydrate former phase and water the various components have different desires to adsorb on liquid water. This depends on the the interaction between each of the hydrate formers and water, as well as the thermodynamic state of the various hydrate formers. In a mixture of CH4 and CO2 then CH4 is superritical and CO2 is subcritical. See for instance Kvamme [22] for an illustration of these aspects; (3) The combined first and second laws of thermodynamics will direct the system toward formation of the most stable hydrates first, under constraints of mass and heat transport. Then there will be a gradual change in hydrate composition.
These aspects are illustrated in Figure 4. Each of these two figures is a combination of two equilibrium calculations and results from dynamic experiments published by Lequang et al. [49]. The experimental data point is obtained using two crystallization methods. High rate means that the experiment is started with high supersaturation. This can be a high ΔP above hydrate stability limit P for a defined temperature. Or it can be large ΔT below hydrate stability T for a fixed pressure. The results from high rate experiment is plotted in Figure 4a together with equilibrium data from this work and from CSMHYD. Low rate experimental data are plotted in Figure 4b with equilibrium calculations. For details on the experiments and timeline for pressure temperature developments the reader is directed to the original study of Lequang et al [49]. In this context the qualitative aspects are the most important.
High thermodynamic driving force witll enhance the impact of points (2) and (3). Criticalpoint for C2H6 is 305.3 K and 48.7 bar and for CO2 the critical point is at 304.1 K and 73.8 bar. The relative driving for for condensation of C2H6 and CO2 versus the supercritical CH4 is substantial. In contrast to the non-polar hydrocarbons the quadropole moment of CO2 also has a favorable interaction with liquid water in the adsorption on water prior to hydrate formation. A fairly flat initial hydrate formation curve indicates dominance of C2H6 and CO2 in the first hydrate films, than the slow crystallization in Figure 4b. Note that the long periods of almost invisible response in the experiments is not nucleation, but very slow transport through hydrate films and thermodynamically controlled rearrangements of hydrate films between hydrate former phase and water. See for instance Kvamme [23] and Kvamme et al. [28] for discussion on nucleation times and critical sizes for CH4 and CO2 hydrates. Critical diameters are on nano scale diameters and associated nucleation times are in the order of nano seconds.
In view of points (1) to (3) and the discussion above hydrates in nature or industry can never reach full thermodynamic equilibrium because there are too many active phases of relevance for hydrate formation and dissociation. Separate fluid phase, liquid water, adsorbed on mineral surfaces and several hydrate phases shows this whether Gibbs phase rule is applied or a full balance of independent thermodynamic variables versus conservations laws or conditions for equilibrium is utilized. As discussed in the next section formation of hydrate from dissolved hydrate formers in water can, mathematically speaking, form infinite number of hydrate phases. All these hydrate will have different composition, different density and different Gibbs free energy. By thermodynamic definition they are therefore all different phases. This can be seen by looking at how the cavity partition function relate to cavity fillings and corresponding mole-fractions in hydrate. In a non-equilibrium system there is no rule that says that chemical potential is the same for all components in all phases. In contrast to an equilibrium system the phase distribution in a non-equilibrium system is determined by point (3) above. Moreover, then also the distribution of each component in each phase related to a unique chemical potential for each components in each phase locally.
ν is fraction of cavity per water. Subscripts large or small means large and small cavity, respectively and i is a guest component index. Corresponding mole-fraction water is then given by: The associated hydrate Gibbs free energy is then:

Hydrate Stability Limits in the Projection of Hydrate Former Concentration in Surrounding Water
Formation of hydrate from solution is possible in between the solubility limit of the actaual guest molecule(s) in water and a lower limit os hydrate stability as function of concentration of the same solutes in water. (18) still applies, but for a defined set of T and P the mole-fraction of hydrate former in the water solution outside the hydrate is now the unknown variable to be solved for in terms of hydrate stability. The actual mole-fraction found in the lower concentration limit for hydrate stability towards water containing hydrate former(s). The relevamt version of (6) is now: in which the superscript aqueous denote chemical potential for the actual hydrate former dissolved in water. For hydrate formers of limited solubility the asymmetric excess convention is the most appropriate to use: The associated ideal gas chemical potential is trivially given by the temperature and the density of the molecule at infinite dilution in water. We have used experimental data for the infinite dilution of methane in liquid water. This is almost constant for variation of pressure and limited dependent on temperature for the relevant conditions. Parameters for the fitted model of activity coefficients are given in Table 4 for Equation (30).  For CO2, a slightly different approach is utilized. The density of CO2 as dissolved in water will correspond to the partial molar volume of CO2 at infinite dilution. The infinite dilution ideal gas chemical potential is not very sensitive to pressure, so the following approximation to only temperature dependency is considered as adequate: where TR is reduced temperature and defined as actual T in Kelvin divided by critical temperature for CO2 (304.35 K). The lower summation 1, 2 indicates starting from 1 and counting in steps of 2. Parameters are given in Table 5 below. The vector sign on mole-fraction x denote mole-fractions and any arrow on top of x denote the vector of all mole-fractions in the actual phase.
Since the chemical potential of CO2 is not necessarily the same for dissolved CO2 in water and CO2 in gas (or liquid) in a non-equilibrium situation, then hydrate formed according to Equation (2) will be different from the first hydrate and accordingly denoted H2. The composition of this hydrate will be different as seen from the corresponding compositions, which follows from Equations (23)- (25).
Some comparisons with experimental data from Yang et al. [50] are shown in Figure 4 for CH4. These comparisons are not directly representative due to the experimental setup and how the experiments are conducted. Furthermore, note that the is a very small pressure dependency in the calculated stability limits, but hardly visible on the scales in Figure 5. These curves are stability limits between two condensed phases. There is a very small poynting correction on the liquid side and a small poynting correction on the hydrate side these will almost cancel. The partial molar volume of water in hydrate is slightly larger than partial molar volume of liquid water, but the impact is not visible over the range of pressures in Figure 5a,b. It is also important to keep in mind that all these calculations are pure predictions. Parameters in the cavity Gibbs free energy of inclusions are derived from Molecular Dynamics simulations and the parameters are the same for all thermodynamic calculations and stability limits in the temperature pressure projection of the stability limit window.
For practical purposes of stability limits towards incoming water through fractures in a sediment and associated dissociation kinetics of in situ hydrate the calculated results are more than accurate enough for the purpose, and the actual error bars related to the experiments are fairly uncertan. Another dimension of this is the hydrate stability window between solubility contour and lowest limit of hydrate stability, which is plotted in Figure 6a. The points for 278.1 K and 278.2 K are hardly visible due to the red contour of the hydrate stability and the yellow contour of the solubility. (a) (b) Figure 6. (a) Calculated pressure temperature contours of minimum liquid mole-fraction CH4 mole fraction for hydrate stability (red) and solubility of CH4 in liquid water (yellow). Experimental data from Yang et al. [50] is plotted as (o)-The values for 278.1 K and 278.2 K are hardly visible here, but can be seen from Figure 6b; (b) mole-fraction CH4 in hydrate formed from saturated CH4 in aqueous solution. Mole-fractions CH4 in hydrate formed along the minimum hydrate stability limits in Figure  4a (red contour) are slightly higher.
We could not find any open hydrate codes based on the difference method that can calculate hydrate formation from various concentrations of dissolved hydrate formers in water. These calculations are critical in many natural situations of hydrates in sediments.
Seawater leaking through fractures will lead to hydrate dissociation due to very limited CH4 in the incoming seawater. This is will lead to fluxes of methane bubbling out from the sediments. Two situations can occur, depending on seafloor temperatures and pressures. If temperature and pressure are outside hydrate stability limits then the methane will bubble out through the seawater column in the fracture. Some CH4 will dissolve and partly be converted through biochemical reactions. If the upcoming gas enters hydrate formation conditions of temperature and pressure at the seafloor, then a dynamically unstable situation develop. Hydrate will dissociate towards seawater under saturated with CH4 from the seafloor side and reform quickly from the reservoir side. This is a complex system in which also the dynamics of biologic comsumption of released CH4 from the top-side can be very important. In some cases the dynamics of geo-biologic ecosystems can even dominate the consumption of released CH4.
In a bigger picture there is a dynamic balance between incoming hydrate formers through fracture systems from below the hydrate formation zones and dissociation of hydrate through fracture systems that brings in seawater from top. If the dissociation flux of hydrate caused by incoming seawater is higher than the flux of new hydrate formation this situation will in the long run lead to depletion of the hydrate. In worst case this can lead to local geo mechanical instabilities and in the worst case landslides.
However, even during hydrate production using for instance pressure reduction hydrate dissociation towards under saturated water can play a significant role. Pumping out water leads to circulation of water from other sections of the sediments through hydrate filled sediments. This incoming water may very well may be water that is under saturated with CH4 and as such phase transitions discussed in this section can assist in hydrate dissociation.
Another important aspects of the residual scheme is that all phases are calculated based on ideal gas as reference state. This results in a very transparent comparison of phase stability which is not possible in the same way with the reference scheme, even if specific parameters are used for also being able to calculate phase transitions discussed in this section. Hydrate compositions and free energies calculated from any route discussed above and below are directly comparable in terms of relative stability. Practically this will be a tool for evaluation of which phases that will dissociate first under various changes of conditions. Moreover, even under constant boundary conditions hydrates of higher Gibbs free energy can be consumed in favor of growth of hydrates of lower Gibbs free energy when supply of new mass is limited.
Gibbs free energy minimization methods for calculating most likely phase distribution and associated compositions, is not trivial within a reference scheme model. The difference in reference level between various phases is one challenge. However, there are also several additional challenges and in the final end it boils down to many parameters that are fitted towards pressure temperature stability limit data.

Thermodynamic Properties
As also discussed above the consistency of free energies for all phases when ideal gas is a universal reference state for all components in all phases is important for calculations of stability limits as well as in kinetic theories on various levels of sophistication, from Classical Nucleation Theory (CNT) [28], Multicomponent Diffuse Interface Theory (MDIT) [51,52], Phase Field Theory (PFT) [53][54][55][56][57][58][59][60][61] or other theories derived from statistical mechanics and concepts from Physics. A brief discussion on free energies are discussed in the next section.
Enthalpy changes related to hydrate phase transitions are needed in any concept for production of CH4 from hydrate. In pressure reduction method the pressure reduction ensures that Gibbs free energy of the system is brought outside of hydrate stability zone, but the enthalpy still has to be primarily supplied from the surrounding formations. Whether the transport capacity and the available heat that can be generated through temperature gradients are sufficient remains to be seen. For the reference method the only possibility is to use the calculated gradients of the pressure temperature in a Clapeyron method as utilized by Anderson [62] or in a simpler scheme as proposed in this work. A much more common approach is the simplified Clausius-Clapeyron which is simplified through neglecting molar volumes of condensed phases. A recently proposed residual thermodynamic scheme for enthalpy calculations [23,[63][64][65][66] is also discussed in the section Enthalpies below. An advantage of the residual thermodynamic scheme for both Free energies and Enthalpies is also the calculation of consistent entropy changes. This is not a key topic in this work, but dynamic entropy generation during various production schemes is an important indicator of production efficiency.

Hydrate Free Energies
An important feature of the residual thermodynamic description for all phases is the possibility to compare stability of various hydrates formed from gas/water, dissolved hydrate formers or adsorbed hydrate former and water. As an example we plot Gibbs free energy contours for CH4 hydrate along the temperature pressure equilibrium curve in Figure 7a below. For comparison we plot free energies of hydrates formed from saturated CH4 in water solution as function of temperature and pressure in Figure 7b. At the saturation limit contours in Figure 7b the chemical potentials for CH4 in the solution is the same as in the gas. Within the limited range of pressures and temperatures the chical potential variation for CH4 is limited and vary around-25 kJ/moles. The similar variation for a three site model of CO2 varies around -39 kJ/moles for the same range of conditions. This difference is the reason for the lower Gibbs free energy of the CO2 hydrates as illustrated in Figure 8a The lower limit of hydrate stability concetration in water surrounding hydrate is found by solving for the concentration of hydrate former in water that gives the same chemical potential for both water and hydrate former in the aqueous phase and the hydrate phase. For CH4 these chemical potentials are close to infinite dilution chemical potential in water and given by Equation (28), with the associated parametrization. For the same ranges of conditions as in Figure 6b the chemical potential of CH4 is in the order of −42.6 kJ/mole. The corresponding free energies for hydrate stability limits in terms of concentration of CH4 varies slightly around −48.6 kJ/mole at 274 K and 50 bar. CO2 is more solvable in water and chemical potential variations larger, but for 274 K and 50 bar the chemical potential of CO2 varies slightly around −48.0 kJ/mole and CO2 hydrate Gibbs free energy varies slightly around −49.3 kJ/mole.

Enthalpies of Hydrate Formation and Dissociation
The Clausius equation for calculating enthalpies of phase transitions is well established and there is no need for a detailed derivation. See for instance Kvamme et al. [66] for a brief review of Anderson's [62] scheme from using the Clapeyron equation: Unlike Anderson's scheme we use Monte Carlo simulation to calculate partial molar volume of guest molecules in the various types of cavities. The Monte Carlo procedures are discussed in much details elsewhere [5,6] and will not be repeated here. The calculated values are listed in Table 6 below and are almost not dependent on temperature for the limited range of hydrate stability in the liquid water region. The molar volume for guest molecules in the gas phase is directly available from the utilized equation of state (SRK). Liquid water molar volume is almost constant and trivially calculated from liquid water density and molecular weight. Hydrate water molecular volume is then calculated according to the following balance for one guest: The density and average molecular weight for hydrate is trivially calculated from the lattice constant (12.01 Å is used as a constant value throughout this work), calculated filling fractions and the corresponding average mole-fractions of water and guest in the hydrate., i.e.: Neglecting volumes of condensed phases volumes (including hydrate volume) reduce Equation (19) to the Clausius-Clapeyron equation: (38) Hydrate formation pressures are generally significant above ideal gas limit for methane and natural gas. Moreover, as discussed above it is not very complicated to calculate the necessary volumes needed for the Clausius equation in (34) so we will not compare results with the Clausius-Clapeyron here.
A fundamentally different approach can be derived from the residual thermodynamics model based on chemical potentials for water and guests [63][64][65][66]: For liquid water, the enthalpy is even more trivially obtained by numerical differentiation of the polyonomial fit of chemical potential as function of T given by Kvamme & Tanaka [4]. For water containinf salt or other compents such as methanol additional contributions according to analyttical or numerical differentiaon of the activity term in Equation (13). An example for methanol addition is given in Figure 9 below.
In an equilibrium situation, the chemical potential of same guest (hydrate former) in the two cavity types must be the same and these must be equal to the chemical potential of the same guest molecule in the phase that it has come from. For the heterogeneous case, this implies chemical potential of the molecule in gas (or liquid) hydrate former phase. However, outside of equilibrium, the gradients in chemical potentials as function of T, P and mole-fractions must reflect how the guest molecule behaves in the cavity.
Enthalpies for various guest molecules in the two types of cavities can be computed by Monte Carlo simulations along the lines described by Kvamme & Lund [5] and Kvamme & Førrisdahl [6] by sampling guest water interaction energies and efficient volumes from the guest molecules movements. That is: U refers to energy and superscript R stands for residual (interaction) contribution. zki denotes compressibility factor for the guest molecule i in cavity k. Consistent ideal gas values for the same interaction models that were applied in evaluation of the residual values is trivial.
where B k means Boltzmann's constant and ki V stands for the excluded volume of a molecule of type i in cavity of type k. This latter volume can be evaluated from the sampled volume of center of mass movements plus the excluded volume due to water/guest occupation. Slightly more complex sampling and calculation for molecules which are not monoatomic (or approximated as monoatomic like methane), but still fairly standard [5,6] and explicit discussion on this is not required here.
For a relevant temperature span in the order of 10 K (273 K-283 K), the differences in enthalpies as evaluated from Equation (40) using Monte Carlo sampled data do not vary substantially and could even be approximated as constant for the purpose of this work. This is as expected because the hydrate water lattice is fairly rigid and the average movements are almost the same for the limited temperature range. Sampled cavity partition functions will of course vary remarkably over the same temperature range because of the direct exponential (Boltzmann factor) dependency. The interaction models for methane (CH4) and carbon dioxide (CO2) used is the same as those used by Kvamme & Tanaka [4]. In addition, note that while there is an average attraction also for carbon dioxide (CO2), the sampled Langmuir constant is very small and not substantial. This is also confirmed by the Molecular Dynamics (MD) studies along the lines of Kvamme & Tanaka [4] whereby the movements of carbon dioxide in the small cavity interferes with several water liberation frequencies and the resulting Gibbs free energy of inclusion is not favorable for carbon dioxide in the small cavity. While small cavity occupation of carbon dioxide has been found at extreme conditions in the ice range of temperatures in some studies [69], it remains unclear if there would be any substantial small cavity filling at all for temperatures above zero degrees Celsius.
The most general approach for calculating enthalpy changes related to temperature pressure stability projection of the phase transition for hydrate formation and dissociation is clearly the residual thermodynamic scheme. Although we have only demonstrated this for pure components here the formalism is totally general for mixtures as well. However, there are not many available studies for mixtures to compare with so it makes sense to start with pure components. Moreover, since CH4 and CO2 are important in the concept for combined CH4 production from hydrate and safe long terms storage of CO2 [23 and references therein] these data are needed by use and likely others. The Clapeyron scheme by Anderson [52] involves fairly many computational steps since it goes through ice. Anderson's [62] scheme is discussed and compared in more detail elsewhere [66]. A much simpler Clapeyron scheme was proposed in this work. Preliminary comparisons for CH4 and CO2 with experimental data as well as the residual scheme are very promising, except for temperatures higher than around 287 K for CH4 and CO2. As expected the Clausius-Clapeyron scheme is inferior and results deviates significantly even for moderate pressures (30 bar) for both CH4 and CO2.
Another aspect that is worthwhile considering is that even for the two very different hydrate formers the enthalpy change as function of temperature is an almost linear function. Practically this means that the specific heat capacity change for the phase transition is almost the same for all the various pressures along the stability curve. The enthalpy change for a fluid as function of pressure is by definition zero for ideal gas, approximately zero up to moderate densities and again almost zero for liquid density. In summary this leaves a window of fairly high density gas at which the pressure dependency is significant. This is also reflected in the difference between the subcritical CO2 in Figure 11c with its higher density than density of CH4 in Figure 10c. Even if the enthalpies in the plots in Figure 10a-d, are not orthonormal since both T and P vary simultaneously in these plots the consequence of the reflection above indicates that it may be feasible to propose the construction two approximate orthonormal functions as: In the most coarse grain approximation a linear approximation of the enthalpy change from the lowest temperature to the highest temperature in Figure 10a,b would practically imply that the pressure dependency in (42) is approximated to zero. In this case calculations of heats of formation and dissociation outside equilibrium curve will be absolutely trivial and the more rigorous approach discussed by Kvamme [51] and Kvamme et al. [53] will not needed for those cases. Some experimental data are even published without information on pressure, as discussed by Kvamme et al [53]. In general there are many limitations on available experimental data for enthalpies of hydrate phase transitions. An example fit is illustrated for CH4 in Figure 11 below. A similar fit for CO2 is given as Figure 12 below.
(a) (b) (c) (d) Figure 10. (a) Calculated enthalpies of hydrate formation along the pressure temperature hydrate stability limit curve for CH4. Solid is calculated using Equation (39), dashed curve is calculated using the Clausius approach in Equation (34)   phase. For the guest molecules the enthalpy for guests in cavities was discussed above and for the guest molecules in fluid (gas, liquid, supercritical) phase the relevant expression for the pure components is separate phase is: The ideal gas contribution will be the same for the component tin the fluid phase and in the cavities so the only differences are in the residual contribution. In summary the pressure sensitivity in the enthalpy is expected to be limited, but not as small as the pressure part of Equation (42) in the fits presented in Figures 11 and 12. A more detailed analysis of the various contributions in Equation (39) would give a different splitting into two approximate orthonormal functions. Nevertheless-the results from (42) may be accurate enough for many practical purposes relative to other uncertainties. Moreover, as also discussed above available experimental data are very incomplete and inconsistent [51,53].
Another limitation of the Clapeyron scheme is the limitation to temperature pressure gradients, which excludes calculation of released heat during formation of hydrate from dissolved hydrate formation in liquid water. A film of hydrate formed from a separate hydrate former phase and water will rapidly create a mass transport barrier. Formation of hydrate from aqueous solution and particularly towards the existing hydrate film, will release heat that will dynamically interact with the mass transport limited growth [72]. Part of the released heat will distribute rapidly through liquid water below, but some of the released heat will dissociate some of the hydrate film. Nano scale (Molecular Dynamics simulations) and meso scale (Phase Field Theory modeling) [53][54][55][56][57][58][59][60][61] may shed more insight into these aspects.
In Figure 13 we plot some calculations of heat of formation for hydrate formed from liquid solution of CO2 in water. Details of the calculation procedures are described in more detail by Kvamme [53]. Basically the calculations follow the same scheme as for hydrate formed from gas hydrate formers and water except that the hydrate former thermodynamics is now a liquid state hydrate former description. The magnitudes of the enthalpies are smaller since the difference between a hydrate former surrounded by more or less structured water in liquid state is closer to the hydrate cavity state. This in contrast to bringing a gas molecule into a cavity of more restricted movements and in most cases higher density. Similar plots for CH4 at the same two temperatures are given in Figure 14 below. The variation in enthalpies for CH4 hydrate formation is smaller for CH4. The different behavior for the two temperatures reflect that the trapping of CH4 in the large cavity of structure I actually involves an expansion for CH4 going from liquid to large cavity. Due to Equations (40) and (41) there is a turning of the gradient in enthalpies of hydrate formation in the pressure region from 100 bar to 200 bar.
While the illustrations for the heterogeneous hydrate formation in Figures 9 and 10 are for hydrate formation along the P, T hydrate stability curve the extension to other conditions is also straightforward as discussed by Kvamme (53) through Taylor expansions. Similar for the homogeneous hydrate formation in Figures 13 and 14.

Discussion
The reference method for calculation of hydrate stability limits in the temperature pressure projection of the hydrate stability window is a very old method with substantial limitations, some of which are illustrated in this work. Present level of molecular simulations, as well molecular models for water and other molecules of relevance for hydrate formation and dissociation, is on a very mature level that should make the reference method redundant. Even the illustrations in this work, using a fairly old water interaction model (the TIP4P potential [73]) illustrates the extended possibilities in thermodynamic analysis related to hydrate phase transitions in natural hydrate systems and industrial hydrates. The fact that natural hydrates in sediments never can reach thermodynamic equilibrium, but reside in a dynamic stationary balance between incoming fluxes of hydrate formers from below and dissociation through fracture systems bringing in water from above that dissociates hydrate should motivate a transition over to the use of residual thermodynamics also for hydrate phases. This will also open up for next generation of hydrate simulators.
The old method of reference properties is substantially weakened by the need for many empirical fittings of fundamental properties (enthalpies, free energies). The focused fit of these properties as well as water/guest interaction fitting limit the use of the reference method to more or less only the pressure temperature projection of the phase stability limits. Other phase transitions, such as dissociation of hydrate towards under saturated water, is not feasible with the old concept. The heat release during hydrate formation and the reverse heat needed to dissociate hydrate, is critical in evaluation of hydrate production. However, enthalpy calculations are also critical in evaluation of dissociation of hydrate plugs and other application related to hydrate phase transitions.
The use of reactive transport analogy and treating each hydrate phase transition as a pseudo reaction gives a totally different platform for hydrate reservoir simulators [74][75][76][77][78] and includes papers in these theses. Hydrate reservoir simulators based on local free energy minimization of competing phase transitions for hydrate formation and hydrate dissociation, under constraints of local mass-and energy-fluxes, has a wider application. A few of these applications have been discussed in this work, but future possibilities include many extensions, as for example couplings to dynamic geo-bio ecosystems.
In addition to a wider application of a residual scheme comes the value of a consistent route to various thermodynamic properties, as illustrated by a new route for calculation of enthalpies of formation and dissociation [51,53,66], as discussed above. This also includes heat of formation from dissolved hydrate formers in water (and corresponding reverse dissociation). We also propose a new and simple Claussius method for environments that do not have a thermodynamic code, but measured pressure temperature stability limit data.
Similar for industrial systems, like processing and transport of hydrocarbon systems, which also has the same situation of not being able to reach equilibrium due to all the possible routes to hydrate formation, including the impact of solid surfaces like rusty pipelines [22,[25][26][27]29,30,[79][80][81]. Even if hydrate forms in a pipeline it can re-dissociate if the flow surrounding the hydrate results in contact with liquid water under saturated with hydrate former or hydrocarbons which is under saturated with water. Moreover, this is of course not limited to hydrocarbons. Any handling of hydrate forming phases that contains water or is flowing together with a water phase in a multiphase pipeline, has to be analyzed in a non-equilibrium fashion.
The reason for the title of the study is a hope that other research groups should start to think about changing from the reference method over to residual thermodynamics. This is also the reason that we provided a very brief discussion of the old method, which is present in many codes around the world today. It is actually very simple to convert codes over to a residual thermodynamic basis. Moreover, there are many reasons for making this change. As we have discussed here some of the advantages of a residual thermodynamic scheme along the lines described here are: (1) The possibility to calculate different hydrate formation and dissociation, phase transitions. This was illustrated for hydrate formation from dissolved hydrate formers and hydrate stability limits; (2) In a general non-equilibrium situation, the advantage of a residual thermodynamic scheme is that every component in every co-existing phase has the same reference state (ideal gas). Direct comparisons of chemical potentials and Gibbs free energies for different phases will therefore also provide a direct comparison of relative phase stabilities and thermodynamic driving forces for phase transitions; (3) Residual thermodynamics link directly into Molecular Dynamics simulations for providing model molecule properties for active phases for which experimental data are impossible to measure. One example is hydrate formers adsorbed on mineral surfaces and subsequent hydrate nucleation toward mineral surfaces. It is possible to measure structures of fluids adsorbed on solid surfaces, but there is no direct coupling over to thermodynamic properties; (4) As illustrated here the advantage is that residual thermodynamic description along the lines described here gives direct and consistent routes to many important thermodynamic properties, as demonstrated with enthalpy of hydrate formation. To our knowledge it is the only available method for calculation of enthalpies of hydrate formation for mixtures; (5) Hydrate nucleation theories are implicit couplings between thermodynamics of the phase Transition (Gibbs free energy change), mass transport dynamics and heat transport dynamics.
All the thermodynamic properties involved in various nucleation theories are available from the concept demonstrated here; (6) Present stage of modeling hydrate production was very limited by lack of consistent thermodynamic tools that is able to address the variety of calculations needed for all the phase transitions involved. The non-equilibrium nature of hydrates in sediments [79][80][81] requires a residual thermodynamic scheme that is able to address competing phase transitions for hydrate formation and hydrate dissociation. Work is therefore in progress [76][77][78][79][80] on the development of a new hydrate reservoir simulator, which is fundamentally. Different from any other hydrate reservoir simulators because in utilize a reactive transport platform in which all hydrate phase transitions are treated as pseudo reactions. Each of the thesis in references [76][77][78][79][80] contain 6 to 12 Journal publications. The thesis can be downloaded from University of Bergen for free or they can be sent from the leading author of this study; (7) The residual thermodynamic scheme described in this work was applied to discussion on maximum water that can be tolerated in various hydrate forming systems during transport in pipelines [22,[25][26][27][28][29][30]. This also includes impact of mineral surfaces (rust) on concentration limit of water in gas before drop-out.
With reference to the title of the study, we have demonstrated that a residual thermodynamic scheme can be a platform for complete thermodynamic description of hydrates in sediments, as well as hydrates forming in industrial situation-To date, we have illustrated this through various hydrate stability limit region, like temperature, pressure and hydrate stability towards surrounding water or gas. Moreover, we have illustrated that the same model can be used to calculate thermodynamic responses like free energy changes and enthalpy changes that are needed in hydrate production and in many other applications. The residual thermodynamic scheme for enthalpies discussed here is quite unique because it can be used for the same multicomponent mixtures as used in other calculations. It is also a consistent scheme since it is derived from the free energy model. Presently the residual scheme as discussed here is the most extensive and general thermodynamic model for hydrate. That does not imply that other researchers need to follow our basic models for residual properties of ice, liquid water and empty hydrate. Molecular Dynamics simulations are very easy today in terms of modern computers, many new models for water-water interactions and many open software packages for conducting the simulations. The message is simply that we should turn over from a limited concept from 1970′s to a more complete thermodynamic model system for more general use. This will open up for a totally different platform that can address many natural systems in a different and more accurate way. This includes dynamic hydrate systems that forms from upcoming gas and dissociates towards incoming seawater and it also includes conventional hydrate seeps that enters seafloor at hydrate forming conditions. The need to understand these systems from a more fundamental thermodynamic point of view is very important in the discussion on changes of carbon fluxes into the oceans. Moreover, hydrate energy is becoming increasingly important for many countries. The residual thermodynamic concept presented here can provide all necessary thermodynamic calculations involved.

Conclusions
The reference method for calculating hydrate stability limits in the temperature pressure projection has many limitations. It is not theoretically sound to fit chemical potentials and enthalpies to a range of different hydrate using a semi-empirical model for the cavity partition functions. While the theoretical platform is a Langmuir type of adsorption theory the semi-empirical aspect comes in how the water lattice is treated as not being disturbed by the guest molecule movements in the cavities. This is fair for small guest molecules like CH4, but may be wrong by one kilojoule per mole for a guest molecule like CO2 in large cavity of structure I. Other semi-empirical aspects are related to various approximations in the evaluation of the Langmuir-constant and the models for interactions between guest molecules and the water molecules in the lattice. Practically all the fitting of parameters in the reference method to a two-dimensional (temperature pressure) projection of the hydrate stability limits is a limitation which makes the concept less useful to address modern hydrate challenges. The rapidly growing interest in hydrate energy requires more accurate thermodynamic description of all dimensions of hydrate stability limits. This involves all phases that can contribute to hydrate formation and hydrate dissociation like aqueous phases and dissolved hydrate formers and adsorption on mineral surfaces. With the rapid development of interaction models for water and other relevant components for hydrate, including mineral surfaces it is time to make more use of molecular dynamics simulations to establish residual thermodynamic models for all phases of relevance for hydrate formation. In this work we have demonstrated that residual thermodynamic modeling for all phases is able to describe a wider range of the hydrate stability limits. Moreover, in addition we have demonstrated that also enthalpies of hydrate formation and dissociation can be predicted by residual thermodynamics. Being able to predict stability limits (free energy related) as well as enthalpies is a good sign of consistency also for entropy development. We have also proposed a promising simple Clapeyron scheme as alternative to other more complex schemes.
The residual thermodynamic scheme presented and illustrated here is totally superior to the old reference method. One of the reasons is that the residual scheme because it provides a consistent scheme for a very wide range of properties that are need in practical applications in natural hydrate systems, as well asand in industrial hydrate systems. This does not mean that other groups need to use our model systems for chemical potentials of water as ice or liquid and water in empty hydrates. The equations that we have presented for hydrate thermodynamic properties, including enthalpy calculations can be applied with any sets of chemical potentials for water derived from molecular modeling.