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

Abstract

In the medium term, gas hydrate reservoirs in the subsea sediment are intended as deposits for carbon dioxide (CO₂) from fossil fuel consumption. This idea is supported by the thermodynamics of CO₂ and methane (CH₄) hydrates and the fact that CO₂ hydrates are more stable than CH₄ hydrates in a certain P-T range. The potential of producing methane by depressurization and/or by injecting CO₂ 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 CO₂ 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 CO₂ is able to build hydrates too, and these hydrates are thermodynamically more stable than methane hydrates. The possibility to destabilize methane hydrate by injecting CO₂ 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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ and the building up and stability of CO₂ hydrates. Uddin et al. [14] used the thermal reservoir simulator STARS to initially simulate the CO₂ 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 CO₂. 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 CO₂ 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=𝐾₀𝑆_𝐺+𝑆_𝐿(1−𝜙)𝑆^𝐺+𝑆_𝐿1−𝜙𝑆_𝐺+𝑆_𝐿^2𝜅𝑘,(3)_rel,𝐺𝑆=𝑘^0rel,𝐺𝑆^∗𝐺,eff^𝜎𝐺,𝑘_rel,𝐿=𝑘^0rel,𝐿𝑆^∗𝐿,eff^𝜎𝐿,(4)_∗𝐺,eff=𝑆^𝐺/𝑆_𝐺+𝑆_𝐿−𝑆^eff_𝐺,res1−𝑆^eff_𝐺,res−𝑆^eff_𝐿,res,𝑆^∗𝐿,_eff=𝑆^𝐿/𝑆_𝐺+𝑆_𝐿−𝑆^eff_𝐿,res1−𝑆^eff_𝐺,res−𝑆^eff_𝐿,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−𝜙𝑆_𝐺+𝑆_𝐿².(6) Here, 𝐾^0 is the permeability under starting conditions including initial hydrate saturation, whereas in (3) 𝐾₀ 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𝑅_CH𝑎_MH𝑐^𝑀,𝐿−𝑐^{∗(𝐻)𝑀,𝐿},=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𝑎_𝛼=⎷𝜙𝑆_𝐺+𝑆_𝐿³2𝐾_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𝑅_CH𝑎^_MH𝑃^𝑀,𝐺−𝑃^{∗(𝐻)𝑀,𝐺},=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−𝜙𝑆_𝐺+𝑆_𝐿^−𝜅𝑆^∗𝐿,eff^−𝜎𝐶.(10) Additionally, the van Genuchten equation is available:𝑃_𝑐=2⋅10⁵𝑆^∗𝐿,eff^−1/0.771−1^1/4.37.(11) The effective saturation 𝑆^∗𝐿,_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) The mass transfer coefficient 𝛽_L is set equal to 2.0⋅10⁻⁸ m/s for all components, gas-liquid interphase area is calculated from𝑎_𝐺𝐿=𝜙𝑆_𝐺6𝑑_𝐵(13) with a constant bubble diameter of 𝑑_𝐵=10⁻⁴ 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.571531+0.003𝜗−1.025⋅10⁻⁵𝜗²𝜂_SW+6.53⋅10⁻¹⁰𝑃_𝐿−0.29𝑠_𝐿,=1.7910−6.144⋅10⁻²𝜗_𝐿+1.4510⋅10⁻³𝜗^2𝐿−1.6826⋅10⁻⁵𝜗^3𝐿−1.5290⋅10⁻⁹𝑃_𝐿+8.3885⋅10⁻¹⁸𝑃^2𝐿+2.4727𝑠_𝐿+𝜗_𝐿+6s.0_L574⋅10⁻¹¹𝑃_𝐿−2.6760⋅10⁻¹⁹𝑃^2𝐿4.8429⋅10⁻²𝜗_𝐿−4.7172⋅10⁻³𝜗^2𝐿+7.5986⋅10⁻⁵𝜗^3𝐿⋅10⁻³.(14) Overall liquid phase heat conductivity and viscosity are calculated with volumetric and molar-weighted functions:ln𝜆𝜂_eff𝐿=(1−𝜙)𝜆_𝑆+𝜙_𝛼≠𝑆𝑆_𝛼𝜆_𝛼,=_𝑖𝑦^𝐿𝑖ln𝜂^𝐿𝑖.(15) Diffusion coefficients are calculated from IFM-GEOMAR data functions [26]:𝛿^SW_𝑀=4.72⋅10⁻¹⁴𝑇_𝐿𝛿_SW𝐶𝜂_𝑊(37.7)^0,6𝜂𝜂^𝑊_SW,=10⁻⁹0.1959+5.089⋅10⁻⁶𝑇_𝐿𝜂_𝑊𝜂𝜂^𝑊_SW,(16) and the Boudreau tortuosity correction [27]:𝛿^eff_𝑖,𝐿=𝛿^𝐿𝑖1−2ln𝜙.(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
𝜕𝜕𝑡𝜙S^𝐺𝜌̃_𝐺+∇⋅𝜌̃_𝐺𝐮_𝐺=1𝑉_𝑖𝜕𝑁_𝑖,𝐺𝜕𝜕𝑡,𝜕𝑡𝜙𝑆^𝐿𝜌̃_𝐿+∇⋅𝜌̃_𝐿𝐮_𝐿=1𝑉_𝑖𝜕𝑁_𝑖,𝐿𝜕𝜕𝑡,𝜕𝑡𝜙𝑆^MH𝜌̃_MH=1𝜕𝑉𝜕𝑁𝜕𝑡^MH,𝜕𝑡𝜙𝑆^CH𝜌̃_CH=1𝑉𝜕𝑁𝜕𝑡^CH.(18)

Mass conservation of components
𝜕𝜕𝑡𝜙𝑆^𝐺𝑐_𝑖,𝐺+∇⋅𝐮_𝐺𝑐_𝑖,𝐺−∇⋅𝜙𝑆_𝐺𝛿_𝑖,𝐺∇𝑐_𝑖,𝐺=1𝜕𝑉𝜕𝑁𝜕𝑡^𝑖,𝐺,𝜕𝑡𝜙𝑆^𝐿𝑐_𝑖,𝐿+∇⋅𝐮_𝐿𝑐_𝑖,𝐿−∇⋅𝜙𝑆_𝐿𝛿_𝑖,𝐿∇𝑐_𝑖,𝐿=1𝑉𝜕𝑁𝜕𝑡^𝑖,𝐿.(19)

Energy conservation for gas, liquid, hydrate, and sediment
𝜕𝜕𝑡𝜙𝑆^𝐺𝜕𝜌̃_𝐺𝑒_𝐺̃+∇⋅𝐮_𝐺𝜌̃_𝐺ℎ_𝐺−∇⋅𝜙𝑆^𝐺𝜆_𝐺∇𝑇_𝐺=𝑞_𝐺̇,𝜕𝑡𝜙𝑆^𝐿𝜕𝜌̃_𝐿𝑒_𝐿̃+∇⋅𝐮_𝐿𝜌̃_𝐿ℎ_𝐿−∇⋅𝜙𝑆^𝐿𝜆_𝐿∇𝑇_𝐿=𝑞_𝐿̇,𝜕𝑡𝜙𝑆^MH−∇⋅𝜙𝜌̃𝑆_MH𝑒_MH̃_MH+𝑆_CH𝜌̃_CH𝑒̃_CH𝜕𝜆_MH+S_CH𝜆_CH∇𝑇_𝐻=𝑞_𝐻̇,𝜕𝑡(1−𝜙)𝜌̃^𝑆𝑒_𝑆̃−∇⋅(1−𝜙)𝜆_𝑆∇𝑇_𝑆=𝑞_𝑆̇.(20)

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−𝑥_𝑖𝑢_𝑅𝑖+1/22⋅𝐵𝑟^𝑖,=𝑢_𝑖+1−𝑢^𝑖+2−𝑢_𝑖+1𝑥_𝑖+2−𝑥_𝑖+1⋅𝑥^𝑖+1−𝑥_𝑖2⋅𝐵𝑟1_𝑖+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 CO₂ 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, CO₂ is injected at the same well for another three years at a maximum rate of 8,000 STD m³/day (Case  2.1).

A variation of the second scenario is the simultaneous production of methane at one well and CO₂ 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 m³/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 CH₄ 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 m³/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 m³ 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 m³ have been recovered, whereas only 50 million STD m³ 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 CO₂ Injection (Case  2.1)

The same reservoir geometry and model parameters as for Case  1 are used to simulate CO₂ 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. CO₂ is injected at 5°C at a maximum rate of 8,000 STD m³/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 CO₂ 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 CO₂ is injected at 5°C into the reservoir, methane hydrate continues decomposing within a larger area. In the region of CO₂ 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 CO₂ 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 CO₂ 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 CO₂ hydrate is strongly limited due to the governing pressure and temperature conditions. As CO₂ migrates into the reservoir, CO₂ 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 CO₂ Injection (Case  2.2)

To simulate production of methane and storage of CO₂ 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, CO₂ at 10°C is injected at a maximum rate of 8,000 STD m³/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 CH₄ 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 CO₂ 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 CO₂ and the heat release due to CO₂ 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 CO₂ injection (approximately 100 m).

After six years, a maximum peak in CO₂ hydrate saturation of about 75% is reached near the CO₂ injection well. Due to unfavourable 𝑃,𝑇 conditions for CO₂ 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 CO₂ 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, CO₂ is sucked towards it. Therefore, more CO₂ 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 CO₂ injection (dashed lines; simulated with the same grid). After 1.5 years, the production rate with CO₂ injection slightly exceeds the results for the case without CO₂ injection. The reason might be that methane is pushed towards the production well while CO₂ is injected. CO₂ 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 m³/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 CO₂ sequestration problem and is set up as a one-well scenario with two consecutive steps (CH₄ hydrate decomposition with methane production/CO₂ hydrate formation) and as a two-well scenario with simultaneous methane production and CO₂ 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, CO₂ is injected with a constant rate of 8,000 STD m³/day in a partially depleted radial reservoir for three years. Within this period, CH₄ hydrate has decomposed further due to unfavourable 𝑃-𝑇 conditions during the injection of CO₂. Simultaneously, CO₂ 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 CH₄ hydrate has decreased to zero, and CO₂ hydrate saturation has increased up to 40% and 75% just next to the well, respectively.

In the case of simultaneous CH₄ production and CO₂ 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 CO₂ leads to an enhanced gas rate (compared to the case without CO₂ injection). A maximum rate of more than 40,000 STD m³/day at the beginning of the production decreases to 35,000 STD m³/day after a period of six years. Within this period, the produced amount of injected CO₂ is almost zero.

The performed simulations show that depressurization with or without simultaneous CO₂ 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.

Symbols

𝑎:	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).

Indices

*:	In equilibrium condition
0:	Intrinsic/entry condition
α:	Phase α
abs:	Absolute
𝐶:	Capillary
𝐶:	Carbon dioxide
CH:	Carbon dioxide hydrate
𝑑:	Diffusive
eff:	Effective
𝐺:	Gas phase
𝐻:	Hydrate phase
𝑖:	Component 𝑖
𝐿:	Liquid phase
𝑀:	Methane
MH:	Methane hydrate
rel:	Relative
res:	Resident
𝑆:	Sediment phase
SW:	Sea water
𝑊:	Pure water
𝑥:	In axial direction
𝑧:	In vertical direction.

Acknowledgments

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.