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

Authors: Georg Janicki; Stefan Schlüter; Torsten Hennig; Hildegard Lyko; Görge Deerberg
Journal: Journal of Geological Research
 (2011)

---

## Abstract

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.


---

## Body

## 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 [2–6]). 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 Figure1 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 (Figure1), 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 Figure2). 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 Table2. 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.


## 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 Figure1 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 (Figure1), 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 Figure2). 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 Table2. 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:(1)uα=-krel⁡,αKabsηα(∇Pα-gρα).
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:(2)krel⁡,α=krel⁡,α(SG,SL),(3)Kabs=K0(SG+SL)((1-ϕ)(SG+SL)1-ϕ(SG+SL))2κ,(4)krel⁡,G=krel⁡,G0(SG,eff*)σG,krel⁡,L=krel⁡,L0(SL,eff*)σL,(5)SG,eff*=(SG/SG+SL)-SG,reseff1-SG,reseff-SL,reseff,SL,eff*=(SL/SG+SL)-SL,reseff1-SG,reseff-SL,reseff.
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:(6)Kabs=K0′(SG+SLSG0+SL0)γ(1-ϕ(SG0+SL0)1-ϕ(SG+SL))2.
Here, K0′ is the permeability under starting conditions including initial hydrate saturation, whereas in (3) K0 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:(7)RMH=1V∂NMH∂t=kMHaMH(cM,L-cM,L*(H)),RCH=1V∂NCH∂t=kCHaCH(cC,L-cC,L*(H)).
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 cM,L*(H),cC,L*(H) 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 kα 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(8)aα=(ϕ(SG+SL))32Kabs(SLSα)2/3,α=(MH,CH).
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:(9)RMH′=1V∂NMH∂t=kMH′aMH′(PM,G-PM,G*(H)),RCH′=1V∂NCH∂t=kCH′aCH′(PC,G-PC,G*(H)).
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 modifiedBrooks-Corey approach also used in [23]:(10)Pc=Pc0((1-ϕ)(SG+SL)1-ϕ(SG+SL))-κ(SL,eff*)-σC.
Additionally, the van Genuchten equation is available:(11)Pc=2⋅105((SL,eff*)-1/0.771-1)1/4.37.
The effective saturation SL,eff* is given by (5). Mass transfer of gaseous components to liquid and vice versa is modelled by a two-film theory linear approach(12)ṅi,GL=βLaGL(ci,G*-ci,G).
The mass transfer coefficient βL is set equal to 2.0·10-8 m/s for all components, gas-liquid interphase area is calculated from(13)aGL=ϕSG6dB
with a constant bubble diameter of dB=10-4 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]:(14)λSW=0.57153(1+0.003ϑ-1.025⋅10-5ϑ2+6.53⋅10-10PL-0.29sL),ηSW=(10-5ϑL31.7910-6.144⋅10-2ϑL+1.4510⋅10-3ϑL2-1.6826⋅10-5ϑL3-1.5290⋅10-9PL+8.3885⋅10-18PL2+2.4727sL+ϑL(6.0574⋅10-11PL-2.6760⋅10-19PL2)+sL(10-5ϑL34.8429⋅10-2ϑL-4.7172⋅10-3ϑL2+7.5986⋅10-5ϑL3))⋅10-3.
Overall liquid phase heat conductivity and viscosity are calculated with volumetric and molar-weighted functions:(15)λeff=(1-ϕ)λS+ϕ∑α≠SSαλα,ln⁡ηL=∑iyiLln⁡ηiL.
Diffusion coefficients are calculated from IFM-GEOMAR data functions [26]:(16)δMSW=4.72⋅10-14TLηW(37.7)0,6(ηWηSW),δCSW=10-9(0.1959+5.089⋅10-6TLηW)(ηWηSW),
and the Boudreau tortuosity correction [27]:(17)δi,Leff=δiL1-2ln⁡ϕ.
Other property data for sediment and hydrates are given in Table 1.Table 1
Property data of hydrates and sediment.
ParameterValueParameterValueHydrate number CH4 hydrate5.90Hydrate number CO2 hydrate7.67Molar mass CH4 hydrate122.3 g/moleMolar mass CO2 hydrate174.3 g/moleDensity CH4 hydrate7,522 mol/m3Density CO2 hydrate6,350 mol/m3Specific heat CH4 hydrate256.8 J/(kg*K)Specific heat CH4 hydrate362.5 J/(kg*K)Thermal conductivity CH4 hydrate0.5 W/(m*K)Thermal conductivity CH4 hydrate0.5 W/(m*K)Dissociation enthalpy CH4 hydrate−54 kJ/moleDissociation enthalpy CO2 hydrate−65 kJ/moleSediment density2,600 kg/m3Sediment thermal conductivity3.9 W/(m*K)Sediment specific heat1,000 J/(kg*K)Table 2
Basic parameters and settings for the depressurization scenario.
ParameterValueParameterValueDiameter of area1,000 mIntrinsic permeabilityK0500 mDHeight of area20 mBrooks-Corey parameterk2,0Porosity of sediment0.50Brooks-Corey parameterσG1.75Initial gas saturation0.05Brooks-Corey parameterσL4.0Initial liquid saturation0.55Effective resident gas saturation0Initial methane hydrate saturation0.40Effective resident liquid saturation0.20Initial pressure77 barsEntry capillary pressure2.8×104 PaInitial temperature10°CCapillary pressure exponentσC0.50Initial salinity (std. sea salt composition)35 g/kg


### 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
(18)∂∂t(ϕSGρ̃G)+∇⋅(ρ̃GuG)=1V∑i∂Ni,G∂t,∂∂t(ϕSLρ̃L)+∇⋅(ρ̃LuL)=1V∑i∂Ni,L∂t,∂∂t(ϕSMHρ̃MH)=1V∂NMH∂t,∂∂t(ϕSCHρ̃CH)=1V∂NCH∂t.Mass conservation of components
(19)∂∂t(ϕSGci,G)+∇⋅(uGci,G)-∇⋅(ϕSGδi,G∇ci,G)=1V∂Ni,G∂t,∂∂t(ϕSLci,L)+∇⋅(uLci,L)-∇⋅(ϕSLδi,L∇ci,L)=1V∂Ni,L∂t.Energy conservation for gas, liquid, hydrate, and sediment
(20)∂∂t(ϕSGρ̃GẽG)+∇⋅(uGρ̃Gh̃G)-∇⋅(ϕSGλG∇TG)=q̇G,∂∂t(ϕSLρ̃LẽL)+∇⋅(uLρ̃Lh̃L)-∇⋅(ϕSLλL∇TL)=q̇L,∂∂t(ϕ(SMHρ̃MHẽMH+SCHρ̃CHẽCH))-∇⋅(ϕ(SMHλMH+SCHλCH)∇TH)=q̇H,∂∂t((1-ϕ)ρ̃SẽS)-∇⋅((1-ϕ)λS∇TS)=q̇S.Fluid velocityuα of phase α is given by Darcy’slaw as stated in Section 3.2. The conservation laws given above are completed by two equations-of-state:(21)ρ̃G=ρ̃G(PG,TG,ci,G),ρ̃L=ρ̃L(PL,TL,ci,L),
a capillary pressure condition,(22)PL=PG-Pc,(23)∂PL∂t=∂PG∂t-∂Pc∂t=∂PG∂t-(∂Pc∂SG)SL∂SG∂t-(∂Pc∂SL)SG∂SL∂t,
respectively, and the phase closing condition:(24)SG+SL+SMH+SCH=1,(25)∂SG∂t+∂SL∂t+∂SMH∂t+∂SCH∂t=0,
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. Ifjc are convective and jd are diffusive fluxes at a boundary, the general form of these BC is(26)jc+jd=0⟶jd=0,
for “no-flow” boundaries, and(27)jc+jd=(jc+jd)out,
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 Section3.5 has the general form(28)A∂u∂t=f(t,x,∂u∂x,1xcx∂∂x(xcxDx∂u∂x),y,∂u∂y,1ycy∂∂y(ycyDy∂u∂y)),
in which u is the vector of solution variables and both the mass matrix A and the generalized diffusion coefficient matrices Dx,Dy may depend on u. The variable c can take values of c = 0, 1, or 2 indicating the coordinate system in the respective spatial direction (c=0: Cartesian, c=1: cylindrical, c=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 (x,y) 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 inx direction are specified. The expressions in y 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 celli and cell (i+1) (indexed (i+1/2)), flux-limited left (L) and right (R) states of the solution variables are defined by(29)ui+1/2L=ui+ui-ui-1xi-xi-1⋅xi+1-xi2⋅B(ri),ui+1/2R=ui+1-ui+2-ui+1xi+2-xi+1⋅xi+1-xi2⋅B(1ri+1),
with(30)B(ri)=ri+|ri|1+|ri|(vanLeer-Limiter),ri=(ui+1-ui)/(xi+1-xi)(ui-ui-1)/(xi-xi-1).
The first-order spatial derivatives are then approximated by(31)∂u∂x|x=xi={ui+1/2L-ui-1/2Lxi+1/2-xi-1/2,if(x-velocity)i≥0,ui+1/2R-ui-1/2Rxi+1/2-xi-1/2,if(x-velocity)i<0.
The (generalized) second-order spatial derivatives are discretized by standard central differences:(32)1xc∂∂x(xcD∂u∂x)|(t,xi)=c+1xi+1/2c+1-xi-1/2c+1⋅[xi+1/2cDi+1/2(ui+1-uixi+1-xi)-xi-1/2cDi-1/2(ui-ui-1xi-xi-1)].
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:(33)Di+1/2=D(xi+1/2,ui+1/2),
or(34)Di+1/2={D(xi+1/2,ui+1/2L),if(x-velocity)i≤0,D(xi+1/2,ui+1/2R),if(x-velocity)i>0.


## 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:(1)uα=-krel⁡,αKabsηα(∇Pα-gρα).
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:(2)krel⁡,α=krel⁡,α(SG,SL),(3)Kabs=K0(SG+SL)((1-ϕ)(SG+SL)1-ϕ(SG+SL))2κ,(4)krel⁡,G=krel⁡,G0(SG,eff*)σG,krel⁡,L=krel⁡,L0(SL,eff*)σL,(5)SG,eff*=(SG/SG+SL)-SG,reseff1-SG,reseff-SL,reseff,SL,eff*=(SL/SG+SL)-SL,reseff1-SG,reseff-SL,reseff.
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:(6)Kabs=K0′(SG+SLSG0+SL0)γ(1-ϕ(SG0+SL0)1-ϕ(SG+SL))2.
Here, K0′ is the permeability under starting conditions including initial hydrate saturation, whereas in (3) K0 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:(7)RMH=1V∂NMH∂t=kMHaMH(cM,L-cM,L*(H)),RCH=1V∂NCH∂t=kCHaCH(cC,L-cC,L*(H)).
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 cM,L*(H),cC,L*(H) 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 kα 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(8)aα=(ϕ(SG+SL))32Kabs(SLSα)2/3,α=(MH,CH).
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:(9)RMH′=1V∂NMH∂t=kMH′aMH′(PM,G-PM,G*(H)),RCH′=1V∂NCH∂t=kCH′aCH′(PC,G-PC,G*(H)).
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 modifiedBrooks-Corey approach also used in [23]:(10)Pc=Pc0((1-ϕ)(SG+SL)1-ϕ(SG+SL))-κ(SL,eff*)-σC.
Additionally, the van Genuchten equation is available:(11)Pc=2⋅105((SL,eff*)-1/0.771-1)1/4.37.
The effective saturation SL,eff* is given by (5). Mass transfer of gaseous components to liquid and vice versa is modelled by a two-film theory linear approach(12)ṅi,GL=βLaGL(ci,G*-ci,G).
The mass transfer coefficient βL is set equal to 2.0·10-8 m/s for all components, gas-liquid interphase area is calculated from(13)aGL=ϕSG6dB
with a constant bubble diameter of dB=10-4 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]:(14)λSW=0.57153(1+0.003ϑ-1.025⋅10-5ϑ2+6.53⋅10-10PL-0.29sL),ηSW=(10-5ϑL31.7910-6.144⋅10-2ϑL+1.4510⋅10-3ϑL2-1.6826⋅10-5ϑL3-1.5290⋅10-9PL+8.3885⋅10-18PL2+2.4727sL+ϑL(6.0574⋅10-11PL-2.6760⋅10-19PL2)+sL(10-5ϑL34.8429⋅10-2ϑL-4.7172⋅10-3ϑL2+7.5986⋅10-5ϑL3))⋅10-3.
Overall liquid phase heat conductivity and viscosity are calculated with volumetric and molar-weighted functions:(15)λeff=(1-ϕ)λS+ϕ∑α≠SSαλα,ln⁡ηL=∑iyiLln⁡ηiL.
Diffusion coefficients are calculated from IFM-GEOMAR data functions [26]:(16)δMSW=4.72⋅10-14TLηW(37.7)0,6(ηWηSW),δCSW=10-9(0.1959+5.089⋅10-6TLηW)(ηWηSW),
and the Boudreau tortuosity correction [27]:(17)δi,Leff=δiL1-2ln⁡ϕ.
Other property data for sediment and hydrates are given in Table 1.Table 1
Property data of hydrates and sediment.
ParameterValueParameterValueHydrate number CH4 hydrate5.90Hydrate number CO2 hydrate7.67Molar mass CH4 hydrate122.3 g/moleMolar mass CO2 hydrate174.3 g/moleDensity CH4 hydrate7,522 mol/m3Density CO2 hydrate6,350 mol/m3Specific heat CH4 hydrate256.8 J/(kg*K)Specific heat CH4 hydrate362.5 J/(kg*K)Thermal conductivity CH4 hydrate0.5 W/(m*K)Thermal conductivity CH4 hydrate0.5 W/(m*K)Dissociation enthalpy CH4 hydrate−54 kJ/moleDissociation enthalpy CO2 hydrate−65 kJ/moleSediment density2,600 kg/m3Sediment thermal conductivity3.9 W/(m*K)Sediment specific heat1,000 J/(kg*K)Table 2
Basic parameters and settings for the depressurization scenario.
ParameterValueParameterValueDiameter of area1,000 mIntrinsic permeabilityK0500 mDHeight of area20 mBrooks-Corey parameterk2,0Porosity of sediment0.50Brooks-Corey parameterσG1.75Initial gas saturation0.05Brooks-Corey parameterσL4.0Initial liquid saturation0.55Effective resident gas saturation0Initial methane hydrate saturation0.40Effective resident liquid saturation0.20Initial pressure77 barsEntry capillary pressure2.8×104 PaInitial temperature10°CCapillary pressure exponentσC0.50Initial salinity (std. sea salt composition)35 g/kg


## 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
(18)∂∂t(ϕSGρ̃G)+∇⋅(ρ̃GuG)=1V∑i∂Ni,G∂t,∂∂t(ϕSLρ̃L)+∇⋅(ρ̃LuL)=1V∑i∂Ni,L∂t,∂∂t(ϕSMHρ̃MH)=1V∂NMH∂t,∂∂t(ϕSCHρ̃CH)=1V∂NCH∂t.Mass conservation of components
(19)∂∂t(ϕSGci,G)+∇⋅(uGci,G)-∇⋅(ϕSGδi,G∇ci,G)=1V∂Ni,G∂t,∂∂t(ϕSLci,L)+∇⋅(uLci,L)-∇⋅(ϕSLδi,L∇ci,L)=1V∂Ni,L∂t.Energy conservation for gas, liquid, hydrate, and sediment
(20)∂∂t(ϕSGρ̃GẽG)+∇⋅(uGρ̃Gh̃G)-∇⋅(ϕSGλG∇TG)=q̇G,∂∂t(ϕSLρ̃LẽL)+∇⋅(uLρ̃Lh̃L)-∇⋅(ϕSLλL∇TL)=q̇L,∂∂t(ϕ(SMHρ̃MHẽMH+SCHρ̃CHẽCH))-∇⋅(ϕ(SMHλMH+SCHλCH)∇TH)=q̇H,∂∂t((1-ϕ)ρ̃SẽS)-∇⋅((1-ϕ)λS∇TS)=q̇S.Fluid velocityuα of phase α is given by Darcy’slaw as stated in Section 3.2. The conservation laws given above are completed by two equations-of-state:(21)ρ̃G=ρ̃G(PG,TG,ci,G),ρ̃L=ρ̃L(PL,TL,ci,L),
a capillary pressure condition,(22)PL=PG-Pc,(23)∂PL∂t=∂PG∂t-∂Pc∂t=∂PG∂t-(∂Pc∂SG)SL∂SG∂t-(∂Pc∂SL)SG∂SL∂t,
respectively, and the phase closing condition:(24)SG+SL+SMH+SCH=1,(25)∂SG∂t+∂SL∂t+∂SMH∂t+∂SCH∂t=0,
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. Ifjc are convective and jd are diffusive fluxes at a boundary, the general form of these BC is(26)jc+jd=0⟶jd=0,
for “no-flow” boundaries, and(27)jc+jd=(jc+jd)out,
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 Section3.5 has the general form(28)A∂u∂t=f(t,x,∂u∂x,1xcx∂∂x(xcxDx∂u∂x),y,∂u∂y,1ycy∂∂y(ycyDy∂u∂y)),
in which u is the vector of solution variables and both the mass matrix A and the generalized diffusion coefficient matrices Dx,Dy may depend on u. The variable c can take values of c = 0, 1, or 2 indicating the coordinate system in the respective spatial direction (c=0: Cartesian, c=1: cylindrical, c=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 (x,y) 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 inx direction are specified. The expressions in y 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 celli and cell (i+1) (indexed (i+1/2)), flux-limited left (L) and right (R) states of the solution variables are defined by(29)ui+1/2L=ui+ui-ui-1xi-xi-1⋅xi+1-xi2⋅B(ri),ui+1/2R=ui+1-ui+2-ui+1xi+2-xi+1⋅xi+1-xi2⋅B(1ri+1),
with(30)B(ri)=ri+|ri|1+|ri|(vanLeer-Limiter),ri=(ui+1-ui)/(xi+1-xi)(ui-ui-1)/(xi-xi-1).
The first-order spatial derivatives are then approximated by(31)∂u∂x|x=xi={ui+1/2L-ui-1/2Lxi+1/2-xi-1/2,if(x-velocity)i≥0,ui+1/2R-ui-1/2Rxi+1/2-xi-1/2,if(x-velocity)i<0.
The (generalized) second-order spatial derivatives are discretized by standard central differences:(32)1xc∂∂x(xcD∂u∂x)|(t,xi)=c+1xi+1/2c+1-xi-1/2c+1⋅[xi+1/2cDi+1/2(ui+1-uixi+1-xi)-xi-1/2cDi-1/2(ui-ui-1xi-xi-1)].
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:(33)Di+1/2=D(xi+1/2,ui+1/2),
or(34)Di+1/2={D(xi+1/2,ui+1/2L),if(x-velocity)i≤0,D(xi+1/2,ui+1/2R),if(x-velocity)i>0.


## 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 Figure6).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.Figures3, 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 t=0, 1, 10, and 36 months.Figure 3
Liquid phase pressurePL history for Case  1.Figure 4
Liquid phase temperatureTL history for Case  1.Figure 5
Methane hydrate saturationSMH 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 Figure6. 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 P,T,s-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. Table3 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).
CaseMass transfer coefficient (m/s)Cumulative production (STD m3)Kinetic constant methane dissociation (m/s)Cumulative production (STD m3)Intrinsic permeability (mD)Cumulative production (STD m3)−2.0E − 93.88E + 072.00E − 104.17E + 072502.22E + 07Normal2.0E − 84.27E + 072.00E − 94.27E + 075004.27E + 07+2.0E − 64.33E + 072.00E − 84.29E + 077506.30E + 07In 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 Table4 are set at the uppermost layer of each zone.Table 4
Basic parameters for simulation of layered reservoir.
Over burdenHydrate layerUnder burdenPorosity of sediment0.50.50.4Intrinsic permeabilityK0′500 mD500 mD500 mDPermeability parameterγ5.85.85.8Initial gas saturation00.050Initial liquid saturation10.551Initial methane hydrate saturation00.400Initial pressure52 bars72 bars77 barsInitial temperature7°C10°C10.6°CInitial salinity35 g/kg35 g/kg35 g/kgThickness100 m20 m100 mAs 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 Figure10). 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 pressurePL history for Case  2.1.Figure 9
TemperatureTL history for Case  2.1.Figure 10
Methane hydrate saturationSMH 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 Figure10). 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 P-T 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.Figure11 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 saturationSCH history for Case  2.1.Figure12 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 thex-axis of the grid in the middle layer in z-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 pressurePL history for Case 2.2.Figure 15
TemperatureTL history for Case 2.2.Figure 16
Methane hydrate saturationSMH history for Case 2.2.Figure 17
Carbon dioxide hydrate saturationSCH history for Case 2.2.The decomposition of methane hydrate takes place at both wells (see Figure16): 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 P,T 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 Figure17 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), Figure19 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.


## 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.Figures3, 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 t=0, 1, 10, and 36 months.Figure 3
Liquid phase pressurePL history for Case  1.Figure 4
Liquid phase temperatureTL history for Case  1.Figure 5
Methane hydrate saturationSMH 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 Figure6. 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 P,T,s-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. Table3 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).
CaseMass transfer coefficient (m/s)Cumulative production (STD m3)Kinetic constant methane dissociation (m/s)Cumulative production (STD m3)Intrinsic permeability (mD)Cumulative production (STD m3)−2.0E − 93.88E + 072.00E − 104.17E + 072502.22E + 07Normal2.0E − 84.27E + 072.00E − 94.27E + 075004.27E + 07+2.0E − 64.33E + 072.00E − 84.29E + 077506.30E + 07In 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 Table4 are set at the uppermost layer of each zone.Table 4
Basic parameters for simulation of layered reservoir.
Over burdenHydrate layerUnder burdenPorosity of sediment0.50.50.4Intrinsic permeabilityK0′500 mD500 mD500 mDPermeability parameterγ5.85.85.8Initial gas saturation00.050Initial liquid saturation10.551Initial methane hydrate saturation00.400Initial pressure52 bars72 bars77 barsInitial temperature7°C10°C10.6°CInitial salinity35 g/kg35 g/kg35 g/kgThickness100 m20 m100 mAs 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 Figure10). 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 pressurePL history for Case  2.1.Figure 9
TemperatureTL history for Case  2.1.Figure 10
Methane hydrate saturationSMH 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 Figure10). 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 P-T 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.Figure11 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 saturationSCH history for Case  2.1.Figure12 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 thex-axis of the grid in the middle layer in z-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 pressurePL history for Case 2.2.Figure 15
TemperatureTL history for Case 2.2.Figure 16
Methane hydrate saturationSMH history for Case 2.2.Figure 17
Carbon dioxide hydrate saturationSCH history for Case 2.2.The decomposition of methane hydrate takes place at both wells (see Figure16): 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 P,T 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 Figure17 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), Figure19 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 Table2), 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 P-T 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.



---

*Source: 462156-2011-10-18.xml*