About this Journal Submit a Manuscript Table of Contents
Journal of Geological Research
Volume 2011 (2011), Article ID 462156, 15 pages
Research Article

Simulation of Methane Recovery from Gas Hydrates Combined with Storing Carbon Dioxide as Hydrates

Fraunhofer Institute for Environmental, Safety, and Energy Technology UMSICHT, Osterfelder Strasse 3, 46047 Oberhausen, Germany

Received 15 June 2011; Accepted 1 August 2011

Academic Editor: Michela Giustiniani

Copyright © 2011 Georg Janicki et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


In the medium term, gas hydrate reservoirs in the subsea sediment are intended as deposits for carbon dioxide (CO2) from fossil fuel consumption. This idea is supported by the thermodynamics of CO2 and methane (CH4) hydrates and the fact that CO2 hydrates are more stable than CH4 hydrates in a certain P-T range. The potential of producing methane by depressurization and/or by injecting CO2 is numerically studied in the frame of the SUGAR project. Simulations are performed with the commercial code STARS from CMG and the newly developed code HyReS (hydrate reservoir simulator) especially designed for hydrate processing in the subsea sediment. HyReS is a nonisothermal multiphase Darcy flow model combined with thermodynamics and rate kinetics suitable for gas hydrate calculations. Two scenarios are considered: the depressurization of an area 1,000 m in diameter and a one/two-well scenario with CO2 injection. Realistic rates for injection and production are estimated, and limitations of these processes are discussed.

1. Introduction

Gas hydrates are ice-like solid compounds of water and gas molecules (clathrates) which are stable at low temperature and elevated pressure [1]. The water molecules build out cages by hydrogen bonds in which gas molecules are embedded. Generally, gas hydrates can contain different guest molecules in different cages, depending on their sizes and the availability of guest molecules under given thermodynamic conditions, but methane is the prevalent gas in natural gas hydrates. The exploitation of natural gas hydrate deposits that are known in various permafrost regions and submarine sediments all over the world is in the focus of several research groups because the amount of methane to be recovered could overcome future energy shortages. The greenhouse gas CO2 is able to build hydrates too, and these hydrates are thermodynamically more stable than methane hydrates. The possibility to destabilize methane hydrate by injecting CO2 as pressurized gas or in liquid form was verified in several small-scale experiments carried out by different research groups (see [26]). Thus, the combination of both processes offers the opportunity to open up new energy resources as well as to combat climate change by reducing CO2 emissions. However, the technical realization of this combination of processes has to face various challenges. Besides the technical and economic efforts for drilling in submarine sediments or in deep layers in permafrost regions, these challenges concern the reaction kinetics and transport resistances within the sediments in which methane hydrates are embedded in natural deposits.

Thus, to find the best strategy for methane recovery from a specific deposit with or without CO2 sequestration, a large variety of parameters describing the properties of the particular layer as well as the time- and location-dependent thermodynamic conditions have to be considered. Within the framework of the German SUGAR (SUbmarine GAs Hydrate Reservoirs) project, strategies to produce natural gas from marine methane hydrates and simultaneously store CO2 as hydrates are explored. Before undertaking drilling tests, numerical simulations of the local processes are necessary and helpful. For this purpose, a new scientific simulation model called UMSICHT HyReS was developed to describe the methane production from submarine hydrate layers and the exchange of methane by carbon dioxide. In addition, the commercially available simulation tool STARS (CMG Ltd., Canada) was used. In the following, the new simulation tool is described, and the results of the calculations based on particular reservoir parameters, reaction kinetics, and extraction techniques are outlined.

2. Principles of Methane Extraction with CO2 Sequestration

2.1. Exploitation Strategies and Their Limitations

All strategies to destabilize methane hydrate are aimed to change the thermodynamic and/or chemical properties within the reservoir in order to overcome the preconditions for stability. The thermodynamic equilibrium curves for methane and carbon dioxide hydrates given in Figure 1 reveal that either thermal stimulation (A-B), depressurization (A-C), or a combination of both approaches may be the best choice to destabilize a hydrate and to recover the pure gas.

Figure 1: Equilibrium curves of methane hydrate and carbon dioxide hydrates for different salinities (stability curves calculated after Tishchenko et al. [7] and with NIST data [8]).

The variety of measures to supply heat into a hydrate layer for thermal destabilization ranges from injecting warm water or steam into a well over injecting fuel and oxygen to especially designed on-site combustion chambers up to electromagnetic heating of a particular region around a well. Recently, Schicks et al. [9] suggested the in situ catalytic oxidation of methane to supply the heat that is needed to decompose the hydrate.

Depressurization of a hydrate-bearing region is achieved by suction through one or more wells. According to the equilibrium curve (Figure 1), lowering the pressure will shift the equilibrium to lower temperatures. However, the release of methane gas from its hydrate consumes heat (e.g., Goel [10]) so that the reaction only continues until the temperature reaches the new equilibrium point. Su et al. [11] reported on depressurization experiments that were carried out in a high-pressure reactor in which methane hydrate had been built within a layer of natural sand. It was shown that the gas production rate does not only depend on the driving force (expressed as the difference of actual and equilibrium pressure). Reaction kinetics, heat and mass transfer limitations, and multiphase flow effects are also important. In addition, the possible formation of ice that may hamper mass and heat transfer must be considered.

2.2. Classification of Hydrate Deposits

In recent years, the classification of natural gas hydrate deposits in three classes according to their geologic and reservoir conditions has become accepted. According to Moridis and Collet [12], class 1 deposits comprise a hydrate layer and an underlying two-phase fluid zone with free gas and formation water. They describe class 2 deposits as a hydrate layer with an underlying mobile water zone with no free gas, and single-hydrate layers without underlying mobile fluids are characterized as class 3 accumulations. Class 1 deposits are the most attractive deposits with respect to exploitation strategies because the bottom of such a layer that is in contact with mobile gas must be in or near thermodynamic equilibrium. Therefore, the energy needed to destabilize the hydrate is minimal compared to classes 2 and 3 in which the hydrate interval may be far within the stability zone. In addition, an underlying layer of mobile water and/or gas provides a better conductivity within the sediment and supplies additional heat. This expands the reaction area where hydrates decompose. Besides these three classes, a fourth class of deposits, that is to say disperse, low-saturated hydrate accumulations in oceanic sediments, was also investigated by Moridis and Sloan [13]. They evaluated an exploitation of such deposits as noneconomic.

Structural properties of geological formations also comprise the conditions for the injection of CO2 and the building up and stability of CO2 hydrates. Uddin et al. [14] used the thermal reservoir simulator STARS to initially simulate the CO2 hydrate formation in four simplified geological formations that were classified with respect to their mean porosity and initial permeability. The parameters (mean porosity 0.25–0.35, permeability 10–1,000 mD) were chosen according to geological formations identified in the Mallik deposit in Canada's Mackenzie Delta which is a permafrost region. The SUGAR project is explicitly aimed at submarine gas hydrates. Hester and Brewer [15] showed over 70 different locations of subaquatic deposits all over the world that had been identified by 2009, with less than half of them being actually known deposits, the others described as inferred deposits that had been identified with indirect markers so far. Some data on the structural properties like mean porosities and intrinsic permeabilities of hydrate-bearing sediments are available for some sites that have been investigated so far (see [14, 16, 17]). Thus mean porosities reported so far range from 0.2 to 0.5 (e.g., [16, 17]), permeabilities range from 0.1 [17] to 1,000 mD [18].

2.3. Simplified Reservoir Model

For simulations, a simplified reservoir model built of layers with homogeneous properties for porosity, permeability, and initial hydrate saturation was used (see Figure 2). According to the classification of Moridis and Collet [12], the reservoir corresponds to class 1 type.

Figure 2: Reservoir model for simulations.

In the simulations presented, the focus is on scenarios where only the hydrate layer with no-flow boundaries for mass and/or heat is considered. The initial conditions are shown in Table 2. As initial pressure, the equilibrium pressure at a given temperature was chosen in order to avoid artificial hydrate formation or decomposition during the simulation run.

3. Simulation Tools

Within the frame of the SUGAR research project, two simulation programs are used:(i)the simulator STARS (Steam, Thermal and Advanced Processes Reservoir Simulator, Computer Modeling Group Ltd. (CMG), Canada), an oil and gas industry standard reservoir simulator,(ii)the simulation code UMSICHT HyReS, a newly developed in-house code.

STARS was primarily developed to model oil and gas flow through 3-dimensional sediments, but it can be set up to deal with hydrate formation and decomposition as well. Due to the complexity of subsea hydrate systems and several restrictions in using STARS regarding full thermodynamics of gas hydrate physics, a second in-house code called UMSICHT HyReS was developed within the project to achieve full flexibility.

3.1. Model Concept

The underlying operation scenario for the development of a new model is methane production from subsea gas hydrate accumulations by depressurization, thermal stimulation, or injection of liquid/soluted CO2. The latter can result in an exchange of methane by carbon dioxide in the hydrate lattice or in a two-step decomposition/formation process with gaseous methane and hydrate-stabilized CO2 in the end [5, 6, 19, 20]. Conservation laws for such a complex situation must be formulated for a gas and a liquid phase flowing through sediment partially blocked by hydrates. Assuming a stable nonmoving sediment phase, only sediment porosity will be constant, whereas the saturation of gas, liquid, and hydrates is free variables. Therefore, the flow within the sediment strongly depends on the time-dependent distribution of phase saturations within the calculated volume, and this is a fundamental feature of the UMSICHT HyReS model. The second fundamental concept is the perception of decomposition/formation kinetics within the hydrates. According to works of Haeckel et al. [21, 22], gas hydrate formation is controlled by a liquid phase formation equilibrium. Conceptually this means that hydrate formation is a two-step process of solution of the gas component in the liquid and formation of hydrate within the liquid. From a thermodynamic point of view, this is more complicated than the well-known kinetics of the Kim-Bishnoi type. The third remarkable feature of the new model is the distinct formulation of energy balances for gas, liquid, hydrate, and sediment phase which was made for dealing with special operation situations in thermally stimulated gas production.

3.2. Flow Models

The velocity of a fluid phase α in the multiphase system is given by Darcy’s law:𝐮𝛼𝑘=rel,𝛼𝐾abs𝜂𝛼𝑃𝛼𝑔𝜌𝛼.(1) Equations of the Carman-Kozeny type are widely used as permeability models for subsea sediments. In UMSICHT HyReS, formulations of Brooks-Corey, van Genuchten, and the model used in STARS are implemented. In this work, the Brooks-Corey model (which is also used by Sun and Mohanty [23]) is applied:𝑘rel,𝛼=𝑘rel,𝛼𝑆𝐺,𝑆𝐿𝐾,(2)abs=𝐾0𝑆𝐺+𝑆𝐿(𝑆1𝜙)𝐺+𝑆𝐿𝑆1𝜙𝐺+𝑆𝐿2𝜅𝑘,(3)rel,𝐺=𝑘0rel,𝐺𝑆𝐺,e𝜎𝐺,𝑘rel,𝐿=𝑘0rel,𝐿𝑆𝐿,e𝜎𝐿𝑆,(4)𝐺,e=𝑆𝐺/𝑆𝐺+𝑆𝐿𝑆e𝐺,res1𝑆e𝐺,res𝑆e𝐿,res,𝑆𝐿,e=𝑆𝐿/𝑆𝐺+𝑆𝐿𝑆e𝐿,res1𝑆e𝐺,res𝑆e𝐿,res.(5) With (3), the blocking of the flow path by solid gas hydrates is included. A similar formulation for the absolute permeability is used in STARS:𝐾abs=𝐾0𝑆𝐺+𝑆𝐿𝑆0𝐺+𝑆0𝐿𝛾𝑆1𝜙0𝐺+𝑆0𝐿𝑆1𝜙𝐺+𝑆𝐿2.(6) Here, 𝐾0 is the permeability under starting conditions including initial hydrate saturation, whereas in (3) 𝐾0 is the permeability under pure conditions without hydrates. Parameterization of these equations is given later in Table 2. The parameters for (6) used in STARS are set to give equal flow as (3) in UMSICHT HyReS. With the concept of relative permeability, the flow-model includes the time-dependent progress of gas, liquid, and hydrate phase saturations and is able to describe subsea sediment flow under changing conditions during methane gas production.

3.3. Kinetic Model

Decomposition and formation of gas hydrates are modelled by a linear driving force approach. For methane and carbon dioxide hydrate formation, the following gross kinetic is assumed:𝑅MH=1𝑉𝜕𝑁MH𝜕𝑡=𝑘MH𝑎MH𝑐𝑀,𝐿𝑐(𝐻)𝑀,𝐿,𝑅CH=1𝑉𝜕𝑁CH𝜕𝑡=𝑘CH𝑎CH𝑐𝐶,𝐿𝑐(𝐻)𝐶,𝐿.(7) As mentioned before, the kinetic model assumes that the formation process is controlled by the liquid phase concentrations of the hydrate-forming gas. Within this approach, the concentrations 𝑐(𝐻)𝑀,𝐿,𝑐(𝐻)𝐶,𝐿 are liquid phase concentrations of the gas components in thermodynamic equilibrium with the hydrates present in the volume. These equilibrium concentrations are different from the solution equilibria between gas and liquid phase (solubility). The hydrate formation from a component present in the gas phase can be modelled as a two-step mechanism of gas solution (step one) and hydrate formation by (7) (step two). Decomposition is modelled by the same equations, but in opposite direction. A formulation for the thermodynamic equilibrium of solution and hydrate formation of methane and carbon dioxide is given by Tishchenko et al. [7]; kinetic constants 𝑘𝛼 are adapted from data given by Haeckel et al. [21, 22], and interphase areas for methane and carbon dioxide hydrate formation/decomposition are calculated according to Sun and Mohanty [23] by𝑎𝛼=𝜙𝑆𝐺+𝑆𝐿32𝐾abs𝑆𝐿𝑆𝛼2/3,𝛼=(MH,CH).(8) Equation (8) ensures that the hydrate decomposition results in a hydrate-free environment; in case of hydrate formation in a hydrate-free environment, a small initial saturation is a prerequisite. Calculations with the STARS code are done with the well-known Kim-Bishnoi kinetic approach:𝑅MH=1𝑉𝜕𝑁MH𝜕𝑡=𝑘MH𝑎MH𝑃𝑀,𝐺𝑃(𝐻)𝑀,𝐺,𝑅CH=1𝑉𝜕𝑁CH𝜕𝑡=𝑘CH𝑎CH𝑃𝐶,𝐺𝑃(𝐻)𝐶,𝐺.(9) Here, the distance to thermodynamic equilibrium given in gas phase partial pressures is used as the driving force for hydrate formation and decomposition. In this approach, liquid phase concentrations of the hydrate forming gas have no influence on hydrate formation.

3.4. Other Submodels and Property Data

Within the scope of UMSICHT HyReS, a large number of submodels for capillary pressure, property data, thermal and diffusional conductivity, and for heat and mass transfer between coexistent phases are used. Most submodels depend on one or more free solution variables (like temperature, pressure, concentration, or saturation), resulting in a strong numerical coupling of balance equations. In the following, the most relevant models are presented.

Capillary pressure is given by a modified Brooks-Corey approach also used in [23]:𝑃𝑐=𝑃0𝑐𝑆(1𝜙)𝐺+𝑆𝐿𝑆1𝜙𝐺+𝑆𝐿𝜅𝑆𝐿,e𝜎𝐶.(10) Additionally, the van Genuchten equation is available:𝑃𝑐=2105𝑆𝐿,e1/0.77111/4.37.(11) The effective saturation 𝑆𝐿,e is given by (5). Mass transfer of gaseous components to liquid and vice versa is modelled by a two-film theory linear approacḣ𝑛𝑖,𝐺𝐿=𝛽𝐿𝑎𝐺𝐿𝑐𝑖,𝐺𝑐𝑖,𝐺.(12) The mass transfer coefficient 𝛽L is set equal to 2.0108 m/s for all components, gas-liquid interphase area is calculated from𝑎𝐺𝐿=𝜙𝑆𝐺6𝑑𝐵(13) with a constant bubble diameter of 𝑑𝐵=104 m. Property data of gas phase components are taken from NIST Standard Reference Database [8] and are molar-weighted; for the gas phase density, the Peng-Robinson equation is used. Liquid phase property data are also taken from NIST and other available sources; special corrections for saline seawater are used for density (UNESCO Standard Reference Equation), thermal conductivity [24], viscosity [25], and diffusion coefficients [26]:𝜆SW=0.571531+0.003𝜗1.025105𝜗2+6.531010𝑃𝐿0.29𝑠𝐿,𝜂SW=1.79106.144102𝜗𝐿+1.4510103𝜗2𝐿1.6826105𝜗3𝐿1.5290109𝑃𝐿+8.38851018𝑃2𝐿+2.4727𝑠𝐿+𝜗𝐿6.05741011𝑃𝐿2.67601019𝑃2𝐿+sL4.8429102𝜗𝐿4.7172103𝜗2𝐿+7.5986105𝜗3𝐿103.(14) Overall liquid phase heat conductivity and viscosity are calculated with volumetric and molar-weighted functions:𝜆e=(1𝜙)𝜆𝑆+𝜙𝛼𝑆𝑆𝛼𝜆𝛼,ln𝜂𝐿=𝑖𝑦𝐿𝑖ln𝜂𝐿𝑖.(15) Diffusion coefficients are calculated from IFM-GEOMAR data functions [26]:𝛿SW𝑀=4.721014𝑇𝐿𝜂𝑊(37.7)0,6𝜂𝑊𝜂SW,𝛿SW𝐶=1090.1959+5.089106𝑇𝐿𝜂𝑊𝜂𝑊𝜂SW,(16) and the Boudreau tortuosity correction [27]:𝛿e𝑖,𝐿=𝛿𝐿𝑖12ln𝜙.(17) Other property data for sediment and hydrates are given in Table 1.

Table 1: Property data of hydrates and sediment.
Table 2: Basic parameters and settings for the depressurization scenario.
3.5. Governing Equations System

Conservation laws are formulated for a two-phase fluid flow through a stable sediment matrix as a heterogeneous continuum model with advection-diffusion type equations for component mass conservation. In UMSICHT HyReS, there are five coexistent phases:(i)gas phase (consists of methane, water, and carbon dioxide),(ii)liquid phase (consists of water, salt, methane, and carbon dioxide),(iii)methane hydrate phase,(iv)carbon dioxide hydrate phase, (v)sediment.

The unsteady conservation conditions to describe the submarine layer flow of gas and liquid lead to a system of partial differential equations. The solution variables of the model are(i)pressure of gas and liquid phase,(ii)saturation of gas, liquid, methane hydrate, and carbon dioxide hydrate,(iii)temperature of gas, liquid, hydrate, and sediment phase,(iv)concentrations of methane, water, and carbon dioxide in the gas phase,(v)concentrations of water, salt, methane, and carbon dioxide in the liquid phase.

Based on these assumptions and settings, the conservation equations used in the UMSICHT HyReS code are given as follows.

Mass conservation of fluid and hydrate phases

Mass conservation of components

Energy conservation for gas, liquid, hydrate, and sediment

Fluid velocity 𝐮𝛼 of phase 𝛼 is given by Darcy’s law as stated in Section 3.2. The conservation laws given above are completed by two equations-of-state: ̃𝜌𝐺=̃𝜌𝐺𝑃𝐺,𝑇𝐺,𝑐𝑖,𝐺,̃𝜌𝐿=̃𝜌𝐿𝑃𝐿,𝑇𝐿,𝑐𝑖,𝐿,(21) a capillary pressure condition,𝑃𝐿=𝑃𝐺𝑃𝑐,(22)𝜕𝑃𝐿=𝜕𝑡𝜕𝑃𝐺𝜕𝑡𝜕𝑃𝑐=𝜕𝑡𝜕𝑃𝐺𝜕𝑡𝜕𝑃𝑐𝜕𝑆𝐺𝑆𝐿𝜕𝑆𝐺𝜕𝑡𝜕𝑃𝑐𝜕𝑆𝐿𝑆𝐺𝜕𝑆𝐿𝜕𝑡,(23) respectively, and the phase closing condition: 𝑆𝐺+𝑆𝐿+𝑆MH+𝑆CH=1,(24)𝜕𝑆𝐺+𝜕𝑡𝜕𝑆𝐿+𝜕𝑡𝜕𝑆MH+𝜕𝑡𝜕𝑆CH𝜕𝑡=0,(25) respectively.

3.6. Boundary Conditions

Boundary conditions (BCs) are necessary to solve the PDE system given above under specific physical constraints. Generally, two types of BC are used in UMSICHT HyReS,(i)“no-flow” conditions at boundaries impermeable for mass and/or heat flow,(ii)“flow” conditions at boundaries permeable for mass and/or heat flow.

Both BC can be mixed, for example, if a boundary is permeable for heat, but impermeable for mass. If 𝑗𝑐 are convective and 𝑗𝑑 are diffusive fluxes at a boundary, the general form of these BC is𝑗𝑐+𝑗𝑑=0𝑗𝑑=0,(26) for “no-flow” boundaries, and𝑗𝑐+𝑗𝑑=𝑗𝑐+𝑗𝑑out,(27) for “flow” boundaries where flow conditions at the point behind the boundary are indexed by “out.” With the general equations (26) and (27), boundary conditions for the solution vector and its derivatives can be set up for a special physical scenario. In UMSICHT HyReS, all boundaries are outer boundaries of the calculated volume and there exist no inner BC.

3.7. Numerical Solution

The system of partial differential equations given in Section 3.5 has the general form𝐀𝜕𝐮𝜕𝑡=𝑓𝑡,𝑥,𝜕𝐮,1𝜕𝑥𝑥𝑐𝑥𝜕𝑥𝜕𝑥𝑐𝑥𝐃𝐱𝜕𝐮𝜕𝑥,𝑦,𝜕𝐮,1𝜕𝑦𝑦𝑐𝑦𝜕𝑦𝜕𝑦𝑐𝑦𝐃𝐲𝜕𝐮,𝜕𝑦(28) in which 𝐮 is the vector of solution variables and both the mass matrix 𝐀 and the generalized diffusion coefficient matrices 𝐃𝐱,𝐃𝐲 may depend on 𝐮. The variable 𝑐 can take values of 𝑐 = 0, 1, or 2 indicating the coordinate system in the respective spatial direction (𝑐=0: Cartesian, 𝑐=1: cylindrical, 𝑐=2: spherical). Equation (28) exhibits the typical structure of a system of convection-diffusion equations. In UMSICHT-HyReS, the well-established (vertical) method of lines is applied for solving: spatial derivatives are approximated by finite differences on a structured rectangular grid in the (𝑥,𝑦) plane, whereas temporal derivatives are left in their continuous form. The result of this semidiscretization is a large system of ordinary differential equations for the values of the solution variables in each grid point. Highly efficient solvers are available for such systems with banded Jacobian matrix, for example, the differential-algebraic integrator RADAU [28].

Below, the discretization schemes used for the first and (generalized) second-order derivatives in 𝑥 direction are specified. The expressions in 𝑦 direction are treated analogously. As will be discussed later in Section 4, the transport of the conserved quantities is strongly convection-dominated. For this reason, special care must be taken in the choice of the discretization scheme for the first-order spatial derivatives in order to avoid unphysical oscillations in the solution variables. In UMSICHT HyReS, a second-order upwind scheme combined with a flux-limiter method as described in [29] is used.

On the interface between cell 𝑖 and cell (𝑖+1) (indexed (𝑖+1/2)), flux-limited left (𝐿) and right (𝑅) states of the solution variables are defined by𝑢𝐿𝑖+1/2=𝑢𝑖+𝑢𝑖𝑢𝑖1𝑥𝑖𝑥𝑖1𝑥𝑖+1𝑥𝑖2𝑟𝐵𝑖,𝑢𝑅𝑖+1/2=𝑢𝑖+1𝑢𝑖+2𝑢𝑖+1𝑥𝑖+2𝑥𝑖+1𝑥𝑖+1𝑥𝑖21𝐵𝑟𝑖+1,(29) with𝐵𝑟𝑖=𝑟𝑖+||𝑟𝑖||||𝑟1+𝑖||𝑟(vanLeer-Limiter),𝑖=𝑢𝑖+1𝑢𝑖/𝑥𝑖+1𝑥𝑖𝑢𝑖𝑢𝑖1/𝑥𝑖𝑥𝑖1.(30) The first-order spatial derivatives are then approximated by𝜕𝑢|||𝜕𝑥𝑥=𝑥𝑖=𝑢𝐿𝑖+1/2𝑢𝐿𝑖1/2𝑥𝑖+1/2𝑥𝑖1/2,if(𝑥velocity)𝑖𝑢0,𝑅𝑖+1/2𝑢𝑅𝑖1/2𝑥𝑖+1/2𝑥𝑖1/2,if(𝑥velocity)𝑖<0.(31) The (generalized) second-order spatial derivatives are discretized by standard central differences:1𝑥𝑐𝜕𝑥𝜕𝑥𝑐𝐷𝜕𝑢|||𝜕𝑥(𝑡,𝑥𝑖)=𝑐+1𝑥𝑐+1𝑖+1/2𝑥𝑐+1𝑖1/2𝑥𝑐𝑖+1/2𝐷𝑖+1/2𝑢𝑖+1𝑢𝑖𝑥𝑖+1𝑥𝑖𝑥𝑐𝑖1/2𝐷𝑖1/2𝑢𝑖𝑢𝑖1𝑥𝑖𝑥𝑖1.(32) For the calculation of the generalized diffusion coefficients at the cell interfaces, there is a choice whether the arithmetic mean or the flux-limited left/right state of the respective solution variable is used:𝐷𝑖+1/2𝑥=𝐷𝑖+1/2,𝑢𝑖+1/2,(33) or𝐷𝑖+1/2=𝐷𝑥𝑖+1/2,𝑢𝐿𝑖+1/2,if(𝑥velocity)𝑖𝐷𝑥0,𝑖+1/2,𝑢𝑅𝑖+1/2,if(𝑥velocity)𝑖>0.(34)

4. Operation Scenarios and Simulation Results

As mentioned earlier, depressurization and the supply of heat or a combination of both methods are the most obvious strategies to destabilize methane hydrate and to recover the emerging gas. An alternative means is to change the methane concentration in the hydrate-building phase by supplying a methane-free substance to the reservoir, for example, water or carbon dioxide. If carbon dioxide is supplied and the operating conditions are chosen properly, it may even be possible to substitute methane by carbon dioxide in hydrate form.

After conforming the relevant inputs in STARS to those of UMSICHT HyReS, two scenarios were set up for simulation:(1)first scenario: pure depressurization of a reservoir with closed boundaries,(2)second scenario: depressurization combined with CO2 injection.

In the first scenario, the pressure in the centre of a cylindrical reservoir with a diameter of 1,000 m is reduced from equilibrium conditions by approximately 50 bars, and methane gas is produced for a period of three years (Case  1). STARS and UMSICHT HyReS are used to study the effect of a variation in model parameters on methane production. The results obtained from both simulators are compared and found in close agreement (see Figure 6).

The second scenario is a continuation of the first and was simulated with STARS. After three years of methane recovery at a bottom hole pressure of 30 bars, CO2 is injected at the same well for another three years at a maximum rate of 8,000 STD m3/day (Case  2.1).

A variation of the second scenario is the simultaneous production of methane at one well and CO2 storage at a second well in a Cartesian reservoir (Case  2.2). The distance between the two wells is 500 m, and the maximum rate of injection is 8,000 STD m3/day. Both reservoirs (radial and Cartesian) almost have the same volume so that the same initial amount of methane is present.

4.1. First Scenario: Depressurization of a Closed Reservoir (Case  1)

The reservoir volume for the depressurization scenario is a cylindrical region with a diameter of 1,000 m and a height of 20 m. The well is situated in the centre of the region and has a diameter of 1 m. Starting from CH4 hydrate equilibrium conditions at about 77 bars and 10°C, the liquid pressure at the well is lowered to 30 bars. Methane hydrate is dissolved and a flow is induced in the direction to the well. Main initial values and model parameters are given in Table 2.

Figures 3, 4, and 5 show the history of pressure, temperature, and methane hydrate saturation over the radius of the reservoir. Pressure decreases slowly from 77 bars to 30 bars within the reservoir because of the restricted permeability of the hydrate formation. Temperature starts at 10°C and decreases to 1°C (the equilibrium temperature corresponding to 30 bars); in the same period of time, the methane hydrate saturation decreases due to methane production in the volume. The process will finally stop at a lower temperature in equilibrium conditions. The overall simulation time was three years (corresponds to the lowest curve in all diagrams). The curves in all diagrams correspond to times 𝑡=0, 1, 10, and 36 months.

Figure 3: Liquid phase pressure 𝑃𝐿 history for Case  1.
Figure 4: Liquid phase temperature 𝑇𝐿 history for Case  1.
Figure 5: Methane hydrate saturation 𝑆MH history for Case  1.
Figure 6: Production rate and cumulative production of methane for Case  1, comparison of simulators.

The production curves for methane gas obtained with UMSICHT HyReS (bold lines) and STARS (thin lines) are shown in Figure 6. To get a better impression of the curve progression, the production period was set to 30 years instead of 3 years. The maximum production rate of about 40,000 STD m3/day is quickly reached right at the beginning of the depressurization. For the next five years, methane production only slightly decreases; after five years, nearly 70 million STD m3 have been produced at the well (UMSICHT HyReS). Thereafter, the supply of gas from the outer part of the reservoir slows down because the pressure gradient towards the well becomes increasingly smaller. Thus, although after 30 years of production, the gas saturation in the reservoir is still at a level of about 12%, the lack of the pressure gradient driving force prevents the free methane gas from being recovered.

Bearing in mind that the models used in UMSICHT HyReS and STARS are of high complexity and differ slightly, the production rate and cumulative gas amount after thirty years are in close agreement. The deviation is mainly caused by different kinetic approaches (see (7)–(9)) and various thermodynamic parameters which are 𝑃,𝑇,𝑠-dependent (as in UMSICHT HyReS) instead of constant values (as in STARS).

Taking the above scenario as the base case, model parameters are varied in order to study the strength of their influence on methane gas production. Table 3 shows the parameters that are varied together with the resulting cumulative methane production after three years. The column in the middle represents the parameter settings of the base case. For the parameter study, the values of the base case are lowered (upper row) and increased (lower row), and the respective cumulative methane production is reported for simulations carried out with UMSICHT HyReS. As can be seen, mass transfer and kinetics are fast enough not to hinder gas recovery. The permeability of the reservoir is found to be the bottleneck of the process as it controls the overall performance of methane production.

Table 3: Influence of transport parameters on methane production (after 3 years).

In addition to the variation of parameters, the boundaries of the reservoir have been extended to include the over and under burdens in the simulation. Both zones are set free of hydrate and gas. As the temperature increases with the depth of the reservoir, a vertical thermal gradient of 30°C/km is assumed. The initial pressure in each layer is taken as the equilibrium pressure corresponding to the given temperature. The values for pressure and temperature listed in Table 4 are set at the uppermost layer of each zone.

Table 4: Basic parameters for simulation of layered reservoir.

As mentioned before, due to the depressurization and the endothermic decomposition of hydrates, the reservoir cools down until a new equilibrium is reached. At this state, no more gas is released from the hydrate, and the production rate decreases. The inclusion of the over and under burdens in the simulation leads to an enhanced gas production: after three years, 132 million STD m3 have been recovered, whereas only 50 million STD m3 are extracted for the base case. The reason is the supply of additional heat from the upper and lower layers. Figure 7 shows that more hydrate decomposes in the vicinity of the layer boundaries than in the middle of the hydrate layer. In order to reproduce that effect in the simulations, the grid near the layers has to be refined. Additionally, a thin impermeable layer was added just above the hydrate layer to avoid the rising of free gas from the hydrate to the upper layers.

Figure 7: Effect of over and under burdens on hydrate decomposition.
4.2. Second Scenario: Depressurization Combined with CO2 Injection (Case  2.1)

The same reservoir geometry and model parameters as for Case  1 are used to simulate CO2 storage after methane production at a single well. The results are obtained with the simulator STARS.

Within the three years of methane production at a well pressure of 30 bars, pressure and temperature in the volume decrease. This pressure decline propagates into the reservoir and causes a decomposition of methane hydrate around the well (see Figure 10). After three years, the well is closed for production and starts to function as an injection well for another three years. CO2 is injected at 5°C at a maximum rate of 8,000 STD m3/day. As can be seen in Figures 8 and 9, pressure and temperature again start to increase close to the well due to the injection process and the exothermic formation of CO2 hydrate.

Figure 8: Liquid phase pressure 𝑃𝐿 history for Case  2.1.
Figure 9: Temperature 𝑇𝐿 history for Case  2.1.
Figure 10: Methane hydrate saturation 𝑆MH history for Case  2.1.

After the three years of depressurization, methane hydrate is decomposed completely in a radius of approximately. 10 m around the well (see Figure 10). As CO2 is injected at 5°C into the reservoir, methane hydrate continues decomposing within a larger area. In the region of CO2 hydrate formation, the 𝑃-𝑇 conditions are below the methane hydrate stability curve. The driving force for methane hydrate decomposition is even larger than in the case of simple depressurization. After the injection phase, the volume free of methane hydrate has a radius of approximately 60 m. In addition, while CO2 is injected, the gaseous methane is pushed away from the well into the (closed) reservoir, and secondary methane hydrate formation occurs because of an increase of pressure due to gas compression.

Figure 11 illustrates the development of CO2 hydrate saturation for the period of injection. Within the first few meters, the saturation reaches values of about 75% after three years of injection. At a distance of more than 10 m away from the well, the formation of CO2 hydrate is strongly limited due to the governing pressure and temperature conditions. As CO2 migrates into the reservoir, CO2 hydrate forms up to 40% within a radius of 55 m around the injection well.

Figure 11: Carbon dioxide hydrate saturation 𝑆CH history for Case  2.1.

Figure 12 illustrates the production and injection curves for Case  2.1 after three years.

Figure 12: Gas production/injection rates and cumulated gas for Case  2.1.
4.3. Second Scenario: Depressurization Combined with CO2 Injection (Case  2.2)

To simulate production of methane and storage of CO2 at two different wells, a Cartesian grid with an edge length of 900 m and a height of 20 m is used. The model parameters and initial conditions are the same as before (see Table 2). All simulations were performed with STARS. Starting from methane hydrate equilibrium at 72 bars (salinity is neglected) and 10°C, the bottom hole pressure in the production well is lowered to 30 bars. Simultaneously, CO2 at 10°C is injected at a maximum rate of 8,000 STD m3/day at the injection well which is located 500 m apart from the production well. Both wells have a diameter of 1.0 m. The geometric configuration is shown in Figure 13.

Figure 13: Reservoir geometry in a Cartesian grid (Case 2.2); half size.

The distributions of pressure, temperature, and hydrate saturations along the 𝑥-axis of the grid in the middle layer in 𝑧-direction are illustrated in Figures 14, 15, 16, and 17, respectively, for a period of six years. As far as the production well (left hand side) is concerned, the situation is the same as for Case 1. The pressure decrease (Figure 14) causes a decomposition of CH4 hydrate and therefore a temperature decline (Figure 15). A front propagates into the reservoir, and gas is produced from hydrates as long as the system is not at equilibrium conditions and a pressure gradient exists. The higher pressure at the injection well (right hand side) leads to the exothermic formation of CO2 hydrate and thus to an increase in temperature around the well (see Figure 15).

Figure 14: Liquid phase pressure 𝑃𝐿 history for Case 2.2.
Figure 15: Temperature 𝑇𝐿 history for Case 2.2.
Figure 16: Methane hydrate saturation 𝑆MH history for Case 2.2.
Figure 17: Carbon dioxide hydrate saturation 𝑆CH history for Case 2.2.

The decomposition of methane hydrate takes place at both wells (see Figure 16): at the production well because of the depressurization and at the injection well because of heat supply by the injection of CO2 and the heat release due to CO2 hydrate formation. As discussed for Case  2.1, the methane hydrate saturation declines due to depressurization within a smaller radius (10 to 20 m) than due to CO2 injection (approximately 100 m).

After six years, a maximum peak in CO2 hydrate saturation of about 75% is reached near the CO2 injection well. Due to unfavourable 𝑃,𝑇 conditions for CO2 hydrate formation near the injection, the saturation decreases significantly within a few meters apart from the well. At a distance of about 20 m, hydrate has only been formed up to a saturation of 45%. The radius within which CO2 hydrate is formed is up to 150 m (see Figure 17).

An effect which can be seen in Figure 17 is clarified in Figure 18. Due to the fact that pressure is lowered at the production well, CO2 is sucked towards it. Therefore, more CO2 hydrate forms between the wells than to the right of the injection well.

Figure 18: Hydrate profile after production and injection period of six years: (a) CH4 hydrate; (b) CO2 hydrate.

In addition to the production rate for Case  2.2 (bold lines), Figure 19 shows the results for a depressurization case without CO2 injection (dashed lines; simulated with the same grid). After 1.5 years, the production rate with CO2 injection slightly exceeds the results for the case without CO2 injection. The reason might be that methane is pushed towards the production well while CO2 is injected. CO2 is also sucked towards the production well, but within the simulated period, the amount reaching the production well is almost zero.

Figure 19: Gas rates and cumulated gas for Case  2.2.

5. Conclusions and Outlook

First, a pure depressurization scenario for a reservoir with 1,000 m in diameter was simulated both with STARS and UMSICHT HyReS. For the given reservoir settings (see Table 2), permeability turns out to be the most important parameter for methane gas production, whereas decomposition kinetics and gas-liquid mass transfer are of lower relevance. A production rate of nearly 40,000 STD m3/day can be achieved over a period of 5 years; over 10 years, the production rate decreases to half of its initial value. The reported rates can be enhanced by means of extended simulation boundaries (including over and under burdens). The results predicted by both simulators are in close agreement because the differences of the models for decomposition kinetics and gas-liquid mass transfer are of lower importance in the cases under consideration. As shown above (see Section 3.3-3.4), the enhanced hydrate formation/decomposition model used in UMSICHT HyReS depends on the concentration of the dissolved gas in the liquid phase and therefore strongly on gas-liquid mass transfer effects. Further simulation cases will be carried out with reservoir data which make mass transfer characteristics more relevant, for example, restricted gas-liquid interphase area. In such cases, model differences between STARS and UMSICHT HyReS will become more obvious.

The second simulation case aims at the industrially relevant CO2 sequestration problem and is set up as a one-well scenario with two consecutive steps (CH4 hydrate decomposition with methane production/CO2 hydrate formation) and as a two-well scenario with simultaneous methane production and CO2 sequestration. These simulations were carried out with STARS in a 2-D radial and 3-D Cartesian volume, respectively. In continuation of the pure depressurization scenario, CO2 is injected with a constant rate of 8,000 STD m3/day in a partially depleted radial reservoir for three years. Within this period, CH4 hydrate has decomposed further due to unfavourable 𝑃-𝑇 conditions during the injection of CO2. Simultaneously, CO2 hydrate has formed around the well. At the end of the operation period, substitution of methane hydrate by carbon dioxide hydrate has occurred only within a radius of about 60 m around the well. The saturation of CH4 hydrate has decreased to zero, and CO2 hydrate saturation has increased up to 40% and 75% just next to the well, respectively.

In the case of simultaneous CH4 production and CO2 storage in a two-well setting with a Cartesian grid, the same effects can be observed as mentioned before. Since the productivity depends on the pressure gradient driving force towards the production well, a local increase of pressure due to injection of CO2 leads to an enhanced gas rate (compared to the case without CO2 injection). A maximum rate of more than 40,000 STD m3/day at the beginning of the production decreases to 35,000 STD m3/day after a period of six years. Within this period, the produced amount of injected CO2 is almost zero.

The performed simulations show that depressurization with or without simultaneous CO2 injection is a practicable way to produce methane from subsea hydrate fields with assured rates over many years. Thermal stimulation is an alternative from a thermodynamic point of view, but there is a lack of appropriate technology to heat up subsea sediment in a wider range around a production well. Some promising ideas (as methane-driven point heaters, e.g.) will be investigated in subsequent work.


𝑎:Volumetric interphase area (m2/m3)
𝐵:van Leer limiter (1)
𝑐:Molar concentration (mole/m3)
̃𝑐𝑝:Molar heat capacity (J/(mole K))
𝐷:Generalized diffusion coefficient (specific)
̃𝑒:Molar total energy (J/mole)
𝑔:Gravity constant (m/s2)
:Molar enthalpy (J/mole)
𝑗:Molar flux (mol/(m2*s))
𝑘:Dissociation constant (m/s)
𝑘rel:Relative permeability (1)
𝐾:Absolute/intrinsic permeability (m2)
𝐿:Length of hydrate layer (m)
𝑁:Molar mass (mole)
̇𝑛:Molar flow per unit volume (mole/(m3s))
𝑃:Pressure (Pa)
𝑃𝐶:Capillary pressure (Pa)
̇𝑞:Heat flow per unit volume (W/(m3s))
𝑅:Volumetric molar conversion rate (mole/(m3s))
𝑆:Saturation = fraction of pore volume (m3/m3)
𝑠:Salinity = mass of salt/mass of pure water (kg/kg)
𝑇:Absolute temperature (K)
𝑡:Time (s)
𝑢:Velocity (volume flow per unit area), solution vector (m/s)
𝑉:Volume (m3)
𝑥:Axial coordinate (m)
𝑦:Molar fraction (mol/mol)
𝑧:Vertical coordinate (m)
β:Mass transfer coefficient (m/s)
𝛿:Diffusion coefficient (m2/s)
𝜀:Volume fraction per unit volume (m3/m3)
𝜙:Porosity (1)
𝛾:Permeability model parameter in STARS code (1)
𝜂:Dynamic viscosity (Pa*s)
𝜅:Brooks-Corey permeability model parameter (1)
𝜆:Heat conductivity (W/(m*K))
𝜈:Hydrate number (1)
𝜗:Celsius temperature (deg C)
𝜌:Specific density (kg/m3)
̃𝜌:Molar density (mole/m3)
σ:Brooks-Corey permeability model parameter (1).


*:In equilibrium condition
0:Intrinsic/entry condition
α:Phase α
𝐶:Carbon dioxide
CH:Carbon dioxide hydrate
𝐺:Gas phase
𝐻:Hydrate phase
𝑖:Component 𝑖
𝐿:Liquid phase
MH:Methane hydrate
𝑆:Sediment phase
SW:Sea water
𝑊:Pure water
𝑥:In axial direction
𝑧:In vertical direction.


The support of our research activities by the German Federal Ministry of Economics and Technology and the German Federal Ministry of Education and Research within the framework of the SUGAR project is gratefully acknowledged.


  1. E. D. Sloan, “Clathrate hydrates: the other common solid water phase,” Industrial and Engineering Chemistry Research, vol. 39, no. 9, pp. 3123–3129, 2000. View at Scopus
  2. M. Ota, Y. Abe, M. Watanabe, R. L. Smith, and H. Inomata, “Methane recovery from methane hydrate using pressurized CO2,” Fluid Phase Equilibria, vol. 228-229, pp. 553–559, 2005. View at Publisher · View at Google Scholar · View at Scopus
  3. M. Ota, K. Morohashi, Y. Abe, M. Watanabe, R. Lee Smith, and H. Inomata, “Replacement of CH4 in the hydrate by use of liquid CO2,” Energy Conversion and Management, vol. 46, no. 11-12, pp. 1680–1691, 2005. View at Publisher · View at Google Scholar · View at Scopus
  4. X. Zhou, S. Fan, D. Liang, and J. Du, “Replacement of methane from quartz sand-bearing hydrate with carbon dioxide-in-water emulsion,” Energy and Fuels, vol. 22, no. 3, pp. 1759–1764, 2008. View at Publisher · View at Google Scholar · View at Scopus
  5. G. Ersland, J. Husebø, A. Graue, B. A. Baldwin, J. Howard, and J. Stevens, “Measuring gas hydrate formation and exchange with CO2 in Bentheim sandstone using MRI tomography,” Chemical Engineering Journal, vol. 158, no. 1, pp. 25–31, 2010. View at Publisher · View at Google Scholar · View at Scopus
  6. H. Lee, Y. Seo, Y. T. Seo, I. L. Moudrakovski, and J. A. Ripmeester, “Recovering methane from solid methane hydrate with carbon dioxide,” Angewandte Chemie, vol. 42, no. 41, pp. 5048–5051, 2003. View at Publisher · View at Google Scholar · View at Scopus
  7. P. Tishchenko, C. Hensen, K. Wallmann, and C. S. Wong, “Calculation of the stability and solubility of methane hydrate in seawater,” Chemical Geology, vol. 219, no. 1–4, pp. 37–52, 2005. View at Publisher · View at Google Scholar · View at Scopus
  8. NIST Chemistry Webbook, http://webbook.nist.gov/chemistry/.
  9. J. M. Schicks, E. Spangenberg, R. Giese, B. Steinhauer, J. Klump, and M. Luzi, “New approaches for the production of hydrocarbons from hydrate bearing sediments,” Energies, vol. 4, no. 1, pp. 151–172, 2011. View at Publisher · View at Google Scholar
  10. N. Goel, “In situ methane hydrate dissociation with carbon dioxide sequestration: current knowledge and issues,” Journal of Petroleum Science and Engineering, vol. 51, no. 3-4, pp. 169–184, 2006. View at Publisher · View at Google Scholar · View at Scopus
  11. K. Su, C. Sun, X. Yang, G. Cgen, and S. Fan, “Experimental investigation of methane hydrate decomposition by depressurizing in porous media with 3-Dimension device,” Journal of Natural Gas Chemistry, vol. 19, pp. 210–216, 2010.
  12. G. R. Moridis and T. S. Collet, “Strategies for gas production from hydrate accumulations under various geological and reservoir conditions,” in Proceedings Tough Symposium, 2003.
  13. G. J. Moridis and E. D. Sloan, “Gas production potential of disperse low-saturation hydrate accumulations in oceanic sediments,” Tech. Rep. LBNL-5268, Lawrence Berkeley National Laboratory, Berkely, Calif, USA, 2006.
  14. M. Uddin, D. Coombe, and F. Wright, “Modeling of CO2-hydrate formation in geological reservoirs by injection of 2 Gas,” Journal of Energy Resources Technology, Transactions of the ASME, vol. 130, no. 3, Article ID 032502, pp. 032502-1–032502-11, 2008. View at Publisher · View at Google Scholar · View at Scopus
  15. K. C. Hester and P. G. Brewer, “Clathrate hydrates in nature,” Annual review of marine science, vol. 1, pp. 303–327, 2009.
  16. S. Hustoft, S. Bünz, J. Mienert, and S. Chand, “Gas hydrate reservoir and active methane-venting province in sediments on < 20Ma young oceanic crust in the Fram Strait, offshore NW-Svalbard,” Earth and Planetary Science Letters, vol. 284, no. 1-2, pp. 12–24, 2009. View at Publisher · View at Google Scholar · View at Scopus
  17. Z. Chen, W. Bai, W. Xu, and Z. Jin, “An analysis on stability and deposition zones of natural gas hydrate in Dongsha Region, North of South China Sea,” Journal of Thermodynamics, vol. 2010, Article ID 185639, 6 pages, 2010. View at Publisher · View at Google Scholar
  18. B. J. Anderson, M. Kurihara, M. D. White et al., “Regional long-term production modeling from a single well test, Mount Elbert Gas Hydrate Stratigraphic Test Well, Alaska North Slope,” Marine and Petroleum Geology, vol. 28, no. 2, pp. 493–501, 2011. View at Publisher · View at Google Scholar
  19. B. A. Baldwin, J. Stevens, J. J. Howard et al., “Using magnetic resonance imaging to monitor CH4 hydrate formation and spontaneous conversion of CH4 hydrate to CO2 hydrate in porous media,” Magnetic Resonance Imaging, vol. 27, no. 5, pp. 720–726, 2009. View at Publisher · View at Google Scholar · View at Scopus
  20. J. W. Jung, D. N. Espinoza, and J. C. Santamarina, “Properties and phenomena relevant to CH4-CO2 replacement in hydrate-bearing sediments,” Journal of Geophysical Research B, vol. 115, no. 10, Article ID B10102, 2010. View at Publisher · View at Google Scholar · View at Scopus
  21. M. Haeckel, E. Suess, K. Wallmann, and D. Rickert, “Rising methane gas bubbles form massive hydrate layers at the seafloor,” Geochimica et Cosmochimica Acta, vol. 68, no. 21, pp. 4335–4345, 2004. View at Publisher · View at Google Scholar · View at Scopus
  22. M. Haeckel, A. Reitz, and I. Klauke, “Methane budget of a large gas hydrate province offshore Georgia, Black Sea,” in Proceedings of the 6th International Conference on Gas Hydrates (ICGH '08), Vancouver, British Columbia, Canada, 2008.
  23. X. Sun and K. K. Mohanty, “Kinetic simulation of methane hydrate formation and dissociation in porous media,” Chemical Engineering Science, vol. 61, no. 11, pp. 3476–3495, 2006. View at Publisher · View at Google Scholar · View at Scopus
  24. D. R. Caldwell, “Thermal conductivity of sea water,” Deep-Sea Research and Oceanographic Abstracts, vol. 21, no. 2, pp. 131–137, 1974. View at Scopus
  25. D. J. Kukulka, B. Gebhart, and J. C. Mollendorf, “Thermodynamic and transport properties of pure and saline water,” Advances in Heat Transfer, vol. 18, pp. 325–363, 1987. View at Publisher · View at Google Scholar · View at Scopus
  26. M. Haeckel, K. Wallmann, et al., “Main equations for gas hydrate modeling,” SUGAR Internal Communication, 2010.
  27. B. P. Boudreau, Diagenetic Models and Their Implementation: Modelling Transport and Reaction in Aquatic Sediments, Springer, Berlin, Germany, 1997.
  28. S. V. Pennington and M. Berzins, “New NAG library software for first-order partial differential equations,” ACM Transactions on Mathematical Software, vol. 20, no. 1, pp. 63–99, 1994. View at Publisher · View at Google Scholar · View at Scopus
  29. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Springer Series in Computational Mathematics 14, Springer, 2nd edition, 1996.