Simulations of long term methane hydrate dissociation by pressure reduction using an extended RetrasoCodeBright simulator

Energy Conversion and Management 68 (2013) 313–323

Contents lists available at SciVerse ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Simulations of long term methane hydrate dissociation by pressure reduction using an extended RetrasoCodeBright simulator Ashok Chejara, Bjørn Kvamme ⇑, Mohammad Taghi Vafaei, Khaled Jemai Department of Physics and Technology, University of Bergen, Allegaten 55, 5007 Bergen, Norway

a r t i c l e

i n f o

Article history: Received 14 May 2012 Received in revised form 26 August 2012 Accepted 3 September 2012 Available online 12 October 2012 Keywords: Methane hydrates RetrasoCodeBright Hydrate dissociation Depressurization Gas production

a b s t r a c t Methane hydrates in sediments are generally not in thermodynamic equilibrium. This implies that there may be several competing hydrate phase transitions in naturally existing hydrate reservoirs. In addition to dissociation due to instability imposed by changes in temperature and/or pressure which can bring hydrate outside stability also gradients in chemical potential caused by concentration gradients may lead to dissociation or formation of hydrate. Mineral surfaces bring additional thermodynamic phases of impact for hydrate phase transitions. The limited numbers of reservoir simulators which have incorporated hydrate are normally simplified by considering only pressure and temperature as criteria for hydrate stability region. In cases where kinetic description is used it is normally based on oversimplified models, typically models derived from experiments in pressure, temperature volume controlled laboratory cells. In the case of hydrate production through pressure reduction heat transport might dominate the kinetics and simplified heat transport kinetic models are frequently in use for this purpose. In lack of reliable data from full scale hydrate production the reservoir simulators are the only tools which can be used to evaluate efficiency of different production scenarios. Several research groups have been recently working on this subject. The approaches for inclusion of hydrates as a phase, and corresponding production simulation results from different simulators vary. In this work we have applied a fundamentally different approach, in which we have reworked a reactive transport reservoir simulator, namely RetrasoCodeBright into a hydrate simulator. This has been accomplished by adding hydrates as pseudo-mineral components. This opens up for non equilibrium thermodynamic description since kinetic models for different competing hydrate phase transitions can be included through their respective kinetic models. The main theoretical tool for generating these kinetic models has been phase field theory simulations, with thermodynamic properties derived from molecular modeling. The detailed results from these types of simulations provide information on the relative impact of mass transport, heat transport and thermodynamics of the phase transition, which enable qualified simplifications for implementation into RetrasoCodeBright. Details of the simulator, and numerical algorithms, are discussed in detail and some relevant examples are shown. In particular we applied the reservoir simulator to data from Mount Elbert methane hydrate deposits from North Slope, Alaska. Pressure reduction is used as gas production method. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction 1.1. Background Amount of natural gas present in the form of gas hydrates is enormous. According to some estimates natural gas hydrates are the largest source of natural gas on earth. Its distribution throughout the world and huge quantity of natural gas trapped in hydrates has resulted in a rapidly increasing interest in hydrate energy [1]. Natural gas hydrate (NGH) is an ice like structure in which water enclathrates small non-polar, and slightly polar, molecules. The ⇑ Corresponding author. Tel.: +47 555 800 00; fax: +47 555 833 80. E-mail address: [email protected] (B. Kvamme). 0196-8904/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enconman.2012.09.003

main component of natural gas is methane and most sources have biogenic origin. Hydrates from thermogenic sources are characterized by the addition of significant amounts of ethane and propane. The two dominating structures for NGH are named structure I (sI) and structure II (sII). Among these types sI is the most common in natural gas hydrate occurrences since biogenic sources of methane dominate. Maximum amount of methane that can occur in a sI methane hydrate is fixed by the clathrate geometry at CH4XH2O where X is less than 5.75 since complete filling of all cavities is impossible to achieve. NGH is found mainly in offshore outer continental margin sediments and in Polar Regions commonly associated with permafrost [2]. NGH is not stable in porous media as a consequence of Gibbs phase rule. The system of fluid/hydrate/ mineral is over determined and several phases compete for the

314

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

Nomenclature F n

p DGi H x

l w

r0 Srl Srg

xh1 ri j0

P0

degree of freedom (–) number of components (–) number of phases (–) Gibbs free energy change (kJ) hydrate phase (–) compositions (mol) chemical potential (kJ/mol) water (–) surface tension at temperature in which P0 has been measured (0.072 Nm1) residual saturation (0) residual saturation (0) mass of solute per mass of liquid (–) effective stress (MPa) pore pressure (MPa)

available mass according to progress towards minimum free energy locally. A typical pitfall in analyses of hydrate stability regions is the unconditional use of the pressure, temperature projection of the phase diagram to construction of hydrate stability regions in reservoir, regardless of flow which my dissociate the hydrate in concentration gradients [3].

1.2. Methane production from natural gas hydrate reservoirs: methods and theoretical studies There are several methods for methane production from natural gas hydrate reservoirs. Depressurization method, in which the pressure is reduced to outside the pressure/temperature stability region for hydrate, is attractive due to the limited associated costs of exploitation. Efficiency varies, depending on a number of different factors like for instance the filling of underlying regions (gas, water, soft clay or shale) and heat transport characteristics of the surroundings. It is currently considered as the most feasible process considering expenses and production rate and has been investigated by many research groups through experiments and simulation studies. Thermal stimulation is another method which has been extensively investigated. Injection of steam or hot water brings the hydrate outside stability region and also supports heat for hydrate dissociation. It is considered to be costly due to huge amount of energy waste to the surroundings. A third method is to use chemicals like for instance thermodynamic hydrate inhibitors, such as methanol or brine, to shift the equilibrium curve and dissociate hydrate. This is a costly option and is hardly considered as a standalone production option. Yet another and more recent novel method is injection of CO2 into the methane hydrate reservoirs. CO2-hydrate is more stable than CH4-hydrate over large regions of temperature and pressure. Mixed CO2/CH4 hydrate, in which CH4 remains trapped in small cavities after conversion, is more stable than either one of the pure hydrates for all pressures and temperatures. Direct conversion of CH4 hydrate into a mixed CO2/CH4 hydrate through CO2 injection is feasible but a slow solid state conversion. Parallel to this new CO2-hydrate formation will be formed from the pore water and injected CO2. The released head from this formation will provide the necessary heat to dissociate in situ methane hydrate. Both mechanisms lead to a win–win situation of combined natural gas production method and CO2 sequestration [4]. During more than two decades several research teams worldwide have conducted theoretical studies of methane production from hydrate reservoirs through modeling on different scales,

Vm

u0 k Sl Sg Pe A P0 k0 Sls Sgs Mi

rij dij

molar volume of a mineral (m3/mol) reference porosity (–) shape function for retention curve (0.457 m) liquid saturation (–) gas saturation (–) CH4 hydrate equilibrium pressure (MPa) constant [1] measured P at certain temperature (0.0196 MPa) intrinsic permeability corresponding to u0 (–) maximum saturation for liquid phase [1] maximum saturation for gas phase [1] molecular weight (g/mol) total stress (MPa) Kronecker symbol (dij = 0 if i – j and dij = 1 if i = j) (–)

including development of hydrate reservoir simulators [5–12]. As a result there are both academic and commercial hydrate simulators available although there are still many uncertainties and questionable approximations in the present generation of simulators. It is beyond the scope of this paper to give a detailed review of these simulators as this is already available through an ongoing code comparison study under the NETL program for methane hydrates [13]. The web-page for this code comparison study is continuously updated and has necessary links to publications on results as well as descriptions of the different simulators. So instead of repeating reviews of these different simulators in this paper we put emphasis on describing an alternative and fundamentally different way of handling the hydrate phases and corresponding phase transition dynamics involved in production modeling. In this paper a different approach according to non-equilibrium nature of hydrate phase transitions in the reservoir will be presented and a new reservoir hydrate simulator will be introduced. This simulator is developed on the RetrasoCodeBright (RCB) [14] reactive transport reservoir simulator platform. An important advantage of this simulator is that the hydrate can be treated as pseudo minerals. This opens up for a transparent implementation of non-equilibrium thermodynamics as competing phase transitions will enter the mass- and energy-balances in a similar fashion as competing reactions of mineral precipitations and mineral dissolutions. The module is designed so that it can easily work according to the non-equilibrium thermodynamic package which is being developed in this group. At this stage kinetic models of hydrate formation and dissociation from phase field theory simulations are used to examine the performance of the module through example cases as mentioned above. 1.3. Theory Common to all current hydrate exploitation simulators is the limitations in thermodynamic description and oversimplified kinetic models. Hydrates in a porous media are generally exposed to at least surrounding aqueous phase and mineral surfaces but also sometimes free gas phase. In the simplest case of pure methane as hydrate former there are only two components and three (or four if free gas is present) phases (aqueous, hydrate, adsorbed). This leaves only one degree of freedom according to Gibbs phase rule and equilibrium cannot be established since both pressure and temperature is defined by hydrodynamics and geothermal gradients.

F ¼npþ2

ð1Þ

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

Eq. (1) simply expresses minimum requirements for a system to be able to reach equilibrium between the p possible phases. It does not say anything about whether these numbers of phases can be reached at the defined conditions of independent thermodynamic variables. For the two-phase system of water and methane outside of hydrate region, there will be two degrees of freedom and equilibrium can be reached at given depths with defined local temperature and pressure. If, on the other hand, the system is inside the hydrate stability zone the system is over determined by one degree of freedom. In this system, as well as other systems which are over or under determined in terms of Gibbs phase rule, directions of minimum free energy development will determine the progress of the system. By considering only pressure/temperature projection and assuming equilibrium according to this projection is in best case an over simplification but it is worse than that. A given system may very well be inside hydrate P–T stability region but will dissociate if the concentrations in neighboring phases are under saturated. This implies for instance that the hydrate might sublimate towards gas if the gas is under saturated with respect to water. Or the hydrate can dissociate towards liquid water phase if the water is under saturated with methane. The system will, however, always tend towards the minimum free energy. When the hydrate is inside its pressure–temperature stability region this means that its free energy is lower than that of the aqueous solution [15]. But the hydrate is only unconditionally stable if the gradient of free energy is favorable (direction of lower free energy) for all independent thermodynamic variables. In addition to temperature and pressure this implies concentrations of all components in all coexisting phases. The term under saturated in a non equilibrium system have a reference to the lowest free energy phase. As an example we can consider carbon dioxide dissolved in water at 1 °C and 150 bars. Extrapolations of Henry’s law would estimate solubility in liquid water of 3.3 mol% carbon dioxide. The same condition with hydrate present will result in 1.6% of carbon dioxide in water co existing with hydrate [16,17]. And if there is no free carbon dioxide phase then equilibrium can be established between hydrate and surrounding aqueous solution. Methane hydrates in geological formations will at least be surrounded by an aqueous phase and hydrate saturation vary substantially, depending on local fluid flows that leads to dissociation. With only hydrate and aqueous phase the degrees of freedom is three and only two independent thermodynamic variables are defined so the system is under determined. Supply of methane from below through fractures and faults add an additional phase and the system might reach equilibrium in the absence of mineral surfaces, which complicates the system further in terms of adsorbed phases which also adds into Gibbs phase rule. Biogenic sources dominate natural gas hydrates and methane concentration is very high. Thermogenic sources of natural gas contain significant amounts of ethane and propane. But this does not change the situation generally because hydrate will first form the most stable form (most propane) and then gradually reduce into methane hydrate when all propane and ethane have been consumed into hydrate. The result is a variation of hydrate phases with different free energies and varying barriers for progress towards more uniform hydrate. In a non equilibrium situation chemical potential for methane is different for different phases. This also implies that hydrate created from gas methane and liquid water will be different from hydrate formed from liquid water and dissolved methane in the water phase. This is visible from the statistical mechanical models for water chemical potential in hydrate, in which the hydrate cavity partition functions are the Langmuir constants times guest fugacity in the ‘‘parent’’ phase. Outside of equilibrium and in a general non equilibrium system the chemical potential of guest molecule in hydrate as well as water chemical potential in hydrate can be

315

estimated by expansion from equilibrium conditions using a Taylor expansion to actual conditions. Mineral surfaces have a dual role in the sense that they will normally impose a lowering of water chemical potential compared to pure liquid water. This implies that hydrate normally cannot form very close to the mineral surfaces. On the other hand the mineral surfaces may serve as adsorption sites for hydrate formers and may as such catalyst heterogeneous hydrate formation slightly (1–2 nm) outside of the mineral surfaces. In summary it appears clear that systems of hydrates in porous media is governed by the possible states of minimum free energy of the system and corresponding kinetics of competing ‘‘reactions’’ as function available mass and heat transport kinetics. Fig. 1 shows a schematic illustration of hydrate in porous media. According to this, three phases of gas, aqueous and hydrate are evident. Methane and water are considered as the only components available in the system. The degree of freedom for such a system is one while it can reduce to 1 if two additional phases of CH4 and water adsorbed to both solid grains and hydrate are also considered, with temperature and pressure fixed locally in the reservoir by hydrostatics or hydrodynamics and geothermal gradient. As a first step into understanding of hydrate formation and dissociation processes based on minimum free energy criteria all possible competing phase transitions must be identified and corresponding free energy changes have to be estimated, Eq. (2) is the free energy change for possible hydrate formation situations. p H;i H;i H;i p DGi ¼ þ½xH;i w ðlw  lw Þ þ xCH4 ðlCH4  lCH4 Þ

ð2Þ

In Eq. (2), ‘H’ represents hydrate phase, ‘i’ represents any of the seven phase transition scenarios, ‘p’ represents liquid, gas and adsorbed phases, ‘x’ composition and ‘l’ chemical potential. The chemical potential for water in the hydrate of type i can be estimated using a modified version of the statistical–mechanical model [18]. Hydrate formation can take place in different case scenarios. Gas and liquid water can form hydrate if temperature and pressure conditions are favorable for it. If liquid is supersaturated with methane, it can also form hydrate homogeneously from aqueous solution [16,17]. Absorbed water and methane gas on mineral surface can also form hydrate [19]. Another possibilities is surface reformation, in this case non-uniform hydrate rearranges due to mass limitations on the interface between an aqueous phase and a gas phase. The initial hydrate film between water and a separate hydrate former phase will rapidly grow to a thickness where mass transport through the hydrate film becomes very small (low diffusivity coefficients). In this case the system will appear in a state of quasi-

I. Hydrate II. Solid Grains III. Empty space is filled with water and gas (CH4) Fig. 1. Schematics of hydrate in porous media.

316

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

equilibrium where the gas side of the hydrate can grow from water in gas to a quasi equilibrium on that side. On the other side of the hydrate film the hydrate can grow from dissolved methane and water. Parallel to this there will be rearrangements of hydrates on the surface since the hydrate will be non-uniform. Hydrates initiates at fairly random spots on the interface and thickest hydrate sections will be more stable than thinner ones in terms of lower free energy. With limited available ‘‘free’’ mass the less stable spots will dissociate in favor of growth of the more stable thick regions. The free energy changes for different ways of creating hydrates described above varies since the chemical potentials of different components are not the same in different co existing phases as long as the system is not able to reach equilibrium (Gibbs phase rule). Eq. (3) below represents hydrate dissociation free energy differences, which is of course in mathematical sense same as Eq. (2) with a different sign in front reflecting that the hydrate is initial phase and fluid phases are resulting products after dissociation. p H;i H;i H;i p DGi ¼ ½xH;i w ðlw  lw Þ þ xCH4 ðlCH4  lCH4 Þ

ð3Þ

Driving forces for hydrate dissociation are condition when pressure and temperature are outside hydrate stability zone, or situations where water is under saturated with methane or sublimation when gas is under saturated with water. Hydrate dissociates into gas and liquid water in case first case. In second case the final phase is water and in the third case (sublimation) final phase will be gas. A non-equilibrium logistic system inside a reservoir simulator analyzes the system locally (every element at every time step). All impossible (DG > 0) and unlikely (|DG|
drate dissociation in this paper. These values for methane hydrate dissociation rates were calculated by Baig [21] using phase field theory. Formation of secondary hydrate was neglected in this study by using hypothetical low value for kinetic hydrate formation rate.

2. RetrasoCodeBright simulator RetrasoCodeBright (RCB) has been used in this study and extended with hydrate phase transitions as ‘‘pseudo reactions’’. RCB is capable of realistic modeling of the reaction rates for mineral dissolution and precipitation, at least to the level of available experimental kinetic data. In contrast to some oil and gas simulators the simulator have flow description ranging from diffusion to advection and dispersion and as such is able to handle flow in all regions of the reservoir, including the low permeability regimes of hydrate filled regions. The mathematical equations for the system are highly non-linear and solved numerically. The numerical approach can be viewed as divided into two parts: spatial and temporal discretizations. Finite element method is used for the spatial discretization while finite differences are used for the temporal discretization. The Newton–Raphson method is adopted for the iterative scheme [26]. RCB is a coupling of a reactive transport code Retraso with a multiphase flow and heat code CodeBright. CodeBright contain an implicit algorithm for solution of flow, heat flow and geomechanical model equations [27–30]. The Retraso extension of CodeBright involves an explicit algorithm for updating the geochemistry [26] as schematically shown in the Fig. 2. Since the simulator, and its algorithms, have been discussed very detailed elsewhere ([27–30] and references therein) we will limit the

Independent Variables(Temperature, Fluid Pressure, Deformation)

CodeBright (Flow, Heat and Geomechanics)

Dependent Variables(Fluxes of Fluids, Saturation, Porosity,...)

Newton-Raphson Iteration

Converged

1.4. Kinetic model In general, as discussed above, hydrates in porous media can never reach equilibrium. It can reach possible states of stationary situations of local minimum free energies under constraints of mass- and heat-flow. All possible phase transition will need some level of super saturation in pressure, temperature and concentrations sufficient to provide nucleation and growth of the new phases. Phase field theory [20] is essentially minimization of free energy under constraints of mass- and heat-transport. Systematic studies of different initial systems and different driving forces will enable a detailed mapping of kinetic rates and rate limiting mechanisms, which facilitate development of simple models suitable for reservoir simulator implementation. Extensive research has been going on in the same group on application of phase field theory in prediction of hydrate formation and dissociation kinetics. See Refs. [21–25] for some examples. Constant kinetic dissociation rate 1  109 (mol/m2/s) has been used for the study of methane hy-

Copy Related Variables to Retraso

Retraso (Reactive Transport Module)

Update Flow Properties, Porosity, Permeability, Salinity etc. Based on Reactive Transport

Newton-Raphson Iteration

Converged

Fig. 2. RCB solves the integrated equations sequentially in one time step.

317

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

discussion of the simulator in this paper instead of extending this paper extensively. RCB is capable of handling both saturated and unsaturated flow, heat transport and reactive transport in both liquid and gas. It is a user friendly code for flow, heat, geomechanics and geochemistry calculation. It offers possibility to just computing the chosen unknowns of user’s interest such as hydro-mechanical, hydro-chemical–mechanical, hydro-thermal, hydro-thermal-chemical–mechanical, and thermo-mechanical. It can handle problems in 1D, 2D and 3D. An important advantage of RCB is the implicit evaluation of geomechanical dynamics. According to Fig. 2, flow, heat and geomechanics are solved initially in CodeBright module through Newton–Raphson iteration and then the flow properties are updated according to the effects of reactive transport on porosity and salinity in a separate Newton– Raphson procedure but for the same time step [26–30]. This implies that the consequences of reactions and pseudo reactions (hydrate phase transitions) are evaluated in a time shifted sequence as function of changes in all phases (phase distributions, compositions, densities) which are transported back to the CodeBright module. CodeBright have been extensively applied with different geomechanical models depending on the type of problems. Basic algorithms involved in CodeBright are described in more details by Olivella et al. [29]. It is beyond the scope and limitations of this paper to evaluate different geomechanical models. This is, however, highly relevant and important and will be subject to separate subsequent studies following the project in this paper. In this paper the geomechanical model applied is a thermo-elastic model [29] and as indicated above we do not make attempt to evaluate the relevance of this model as the main focus have been on the implications of non-equilibrium thermodynamics and application of reactive transport method. This makes it possible to study the implications of fast kinetic reactions such as hydrate formation or dissociation more realistically. Hydrate in the reservoir can form from different phases such as liquid phase, gas phase and adsorbed phase on minerals and according to the phase rule there will be always several competing reactions that will prevent the system from reaching equilibrium. For instance hydrate formation rate from CH4(g) + H2O(l) will not be the same as CH4(aq) + H2O(l). Because of this non-equilibrium nature of hydrate in the reservoir, it is very important to consider the kinetics of hydrate formation and dissociation for all possible scenarios [20]. For this purpose, nonequilibrium thermodynamics of hydrate should be employed to determine the kinetic rates of different competing scenarios in each node and each time step according to the temperature, pressure and composition of the system. In this study CH4 hydrate has been added into the simulator as a pseudo-mineral component with a constant kinetic rate for hydrate formation and dissociation so that hydrate formation and dissociation possibility can be studied in the aquifer. Hydrate formation and dissociation can directly be observed through porosity changes in the specific areas of aquifer. Porosity reduction indicates hydrate formation and porosity increase indicates hydrate dissociation. Temperature, pressure and CH4 concentration in all possible phases are three factors, which influence hydrate formation or dissociation. In future, constant kinetic rates for hydrate reactions will be replaced by a thermodynamic code to account for all different competing reactions. This feature of non-equilibrium evaluation of hydrate distinguishes RCB from all other academic and commercial hydrate simulators available today. To illustrate impact of hydrate formation on the reservoir geomechanical properties a simple 2D model case was constructed. The results from this model for some important mechanical, hydraulic properties were assembled and illustrated through a graphical window post-processor GiD [31]. A brief overview of independent variables, constitutive equations and equilibrium restriction are given in Tables 1 and 2 respectively.

Table 1 Equations and independent variables. Equation

Variable name

Equilibrium of stresses Balance of liquid mass Balance of gas mass Balance of internal energy

Displacements Liquid pressure Gas pressure Temperature

Table 2 Constitutive equations and equilibrium restrictions. Constitutive equation

Variable name

Darcy’s law Fick’s law Fourier’s law Retention curve Mechanical constitutive model Phase density Gas law

Liquid and gas advective flux Vapor and gas non-advective flux Conductive heat flux Liquid phase degree of saturation Stress tensor Liquid density Gas density

Equilibrium restrictions Henry’s law Psychometric law

Air dissolved mass fraction Vapor mass fraction

2.1. The independent variables The governing equations for non-isothermal multiphase flow of liquid and gas through porous deformable saline media have been established. Variables and corresponding equations are tabulated in Table 1. 2.2. Constitutive equations and equilibrium restrictions Associated with this formulation there is a set of necessary constitutive and equilibrium laws. Table 2 is a summary of the constitutive laws and equilibrium restrictions. The dependent variables that are computed using each of the laws are also included. 2.3. Reactive transport in Retraso The Retraso part of the code has a built in state of the art geochemical solver and in addition capabilities of treating aqueous complexation (including redox reactions) and adsorption. Different types of nonlinear models are already implemented in the CodeBright part of the code and the structure of the code makes it easy to implement new models derived from theory and/or experiments. The current version of RCB has been extended from ideal gas into handling of CH4 according to the Soave–Redlich–Kwong (SRK) equation of state. This equation of state is used for density calculations as well as the necessary calculations of fugacity of the CH4 phase as needed in the calculation of dissolution of CH4 into the water. The dissolution and precipitation of minerals, with the exception of some carbonates, is slow natural processes and in our example mineral/fluid reactions are described by built-in kinetics based on available experimental data from open literature [21]. 2.4. Permeability and porosity calculations in RCB Reactive transport properties can also affect the thermo hydraulic problem. RCB can model the effect of dissolution and precipitation of minerals on porosity and permeability. The change in porosity (Du) is calculated from the change in the concentrations of the minerals (Dcm) through:

Du ¼ 

X V m;i Dcm;i i

ð4Þ

318

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

Intrinsic permeability and relative permeability of liquid/gas are calculated based on porosity. For calculation of intrinsic permeability, Kozeny’s model is used [32].

k ¼ ko 

u3 ð1  uo Þ2  u3o ð1  uÞ2

ð5Þ

Liquid phase relative permeability is given as

krl ¼ A  ðSe Þk

ð6Þ

Gas phase relative permeability is given as

krg ¼ A  ðSeg Þk

ð7Þ

where Se is effective saturation for liquid and Seg is effective saturation for gas. These are calculated from van Genuchten model [33]:

"   1 #k Sl  Srl Pg  Pl 1k Se ¼ ¼ 1þ Sls  Srl P Seg ¼

ð8Þ

Sg  Srg Sgs  Srg

ð9Þ

r r0

ð10Þ

P ¼ P0 

3.1. CH4 hydrate equilibrium condition (equilibrium Pe vs T relationship) We used CH4 hydrate equilibrium hydration pressure Pe vs temperature (T > 0 °C) relationship proposed by Sloan [36] and Alp et al. [37]. It describes the conditions of coexistence of CH4-hydrate, gas and water, and is given by Eq. (13) below:

LnPe ¼ 1:94138504464560  105 þ 3:31018213397926  103  T

The surface tension (r) is calculated as given in Eq. (11) below [34]:



252:93 r ¼ 0:03059  exp 273:15 þ T

xhl ¼

elements. In total there are 1312 such elements and 1386 nodes in it. The model consisted of three layers; one hydrate layer in the middle and two ‘impermeable’ shale layers on top and below the hydrate layer. Hydrate saturation in hydrate layer was equal to 65% by volume initially. Hypothetical values for mineral compositions and properties in shale layers have been considered. The reason for doing is that shale minerals like quartz and calcites have low dissociation rates compare to methane hydrate. So, dissociation of shale mineral is negligible compare to dissociation of methane hydrate. Thickness of hydrate layer is 12 m and thickness of shale layers are 30 m each. Gas production well is located in the left side of the reservoir and it is 0.2 m wide.



þ 0:04055  x

P i cai M i P 1 þ i cai M i

h l

ð11Þ

 2:25540264493806  101  T 2 þ 7:67559117787059  102  T 3  1:30465829788791  104  T 4 þ 8:86065316687571  108  T 5

ð13Þ

Reservoir properties are given in Tables 3–5. These tables contain data related to chemical species, reservoir properties and material properties.

ð12Þ 4. Long term simulation results and discussion

3. Description of two dimensional model A hypothetical reservoir model based on class 3 hydrate deposits at Mount Elbert, Alaska, has been built. Geometry of the 2D model is shown in Fig. 3. Some reservoir properties such as temperature, pressure, depth, porosity and permeability have been taken from original Mount Elbert hydrate reservoir site C-unit [35]. Temperature and production well pressure have been kept constant. The geometry of the 2D domain is 450 m  72 m rectangle. For simulation purpose this reservoir has been divided into smaller

Results for hydrate saturation, gas production rate, gas pressure and stresses in xx, yy directions have been generated in RCB simulator at regular time intervals. Some of the results are presented and discussed in this paper with the help of graphical representations of outputs. 4.1. Spatial distribution of hydrate saturation Spatial distribution of hydrate saturation is a visualization of hydrate dissociation development in the hydrate layer with time.

Fig. 3. Geometry of 2D model with CH4 production well.

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323 Table 3 Chemical species in different forms. Species Aqueous

Hydrate layer 

H2O, OH ,

O2(g), CO2 3 ,  H2 SiO2 4 ; HSiO3

2 H2O, OH, HCO 3 , O2(g), CO3 ,

H , SiO2(aq), CH4 hydrate (65%)

 H+, SiO2(aq), H2 SiO2 4 ; HSiO3 CH4 hydrate (0%)

Quartz CH4

Quartz –

+

Hydrate saturation Rock mineral Gas

Shale layers

HCO 3,

Table 4 Reservoir properties. Parameter

Hydrate layer

Shale layers

Pressure (MPa) Temperature (°C) Production well pressure (MPa)

7.2 4.0 2.7

7.2 4.0 –

Table 5 Material properties.

319

fluxes have been converted into gas flow rates at standard temperature and pressure conditions. Fig. 5 shows simulated gas production rates with time. As can also be observed from Fig. 5 the initial gas production is small. After few years of low production a significant rise has been observed after roughly 5 years. Since there was no original gas in the reservoir the gradual increase in released gas and corresponding two phases clearly affects the overall gas production. Movement of gas is dependent on porosity and permeability of the reservoir. As discussed above permeability vs hydrate saturation needs to be reworked and this also goes for the relative gas/liquid permeability correlations. As solid hydrate dissociates, porosity and permeability of hydrate layer increases according to Eqs. (4) and (5) respectively. This increase in permeability assists in fluid flow. During total production period of 50 years, maximum volumetric production rate achieved has been close to 8800 (ST m3/day). It is clear that production rates up to 10 years have been low. Increase in production rates has been relatively steep between 10 years and 30 years. The subsequent increase after that time showed a reducing trend. 4.3. Spatial distribution of pressure

Property

Hydrate layer

Shale layers

Young’s modulus (GPa) Poisson’s ratio Zero stress porosity Permeability (mD) (with hydrate) Permeability (mD) (without hydrate) van Genuchten’s gas entry pressure (kPa), (at zero stress) van Genuchten’s exponent (m)

0.5 0.25 0.35 0.12 1000 196

0.5 0.25 0.35 0 0 196

0.457

0.457

Available fluid volume fraction of the porosity (as defined by total available fraction pore volume) is the apparent porosity available for fluid flow. RCB updates porosity according to mineral dissociation or precipitation and this criterion has also been used in same fashion to calculate hydrate saturation based on porosity change. The correlations applied need to be revised and reworked in further development in future. Fig. 4 illustrates spatial distribution of hydrate saturation at different time steps. Pressure of gas production well has been kept constant at 2.7 MPa for the entire simulation. This low pressure resulted in a subsequent pressure drop in hydrate layer and resulting hydrate dissociation. The corresponding reduction of hydrate saturation can be viewed in Fig. 4 for different time steps. Changes in saturation have been slow initial but significant increase has been observed after the first 10 years. Hydrate dissociation is a function of temperature, pressure and concentrations of water and CH4 gas. Initially low temperature conditions in reservoir prevented hydrate from dissociating. This also reflects in the results which show that changes in hydrate saturation are not significant until first 3 years. Only small changes are noticeable in the left side of hydrate layer near production well. It can also be concluded that dissociation started near production well and propagated towards right side as time progressed. Results also indicates that regional disparities in hydrate saturation are present in earlier results from 3 and 10 years as illustrated in Fig. 4a and b. These disparities vanish afterward and homogeneous spatial distribution of hydrate saturation in entire hydrate layer has then been observed after significant time as illustrated in Fig. 4c and d respectively. 4.2. Gas production rates RCB generates output results corresponding to fluid flows for gas and liquid in terms of respective phase fluxes (m/s). Gas phase

Spatial distribution of pressure in hydrate layer is illustrated in Fig. 6. Initial pressure at every node in reservoir was 7.2 MPa. Introduction of a low pressure well had significant impact on reservoir pressure. After 10 years, pressure in hydrate layer got reduced to significantly lower values that generated favorable conditions for hydrate dissociation. Low initial permeability reduced the response on reservoir pressure and a corresponding slow reduction was observed. Eventually this reduction spread over the entire hydrate layer in the period between 1 month and a year. For non-permeable shale layer, very low value of permeability (negligible compare to hydrate layer) has been used. This has protected shale layer from any pressure reduction and negligible fluid flow has been observed in shale layer. 4.4. Stresses The implicit geomechanical module in RCB has the advantage of no time shift between flow analysis and geomechanical impact. The geomechanical analysis does not involve any trigonometric functions and the effective stress in each element can be directly compared to the tensile strength of the material. Study of stresses in reservoir is important to predict how reservoir is going to respond to the volume and density changes involved in gas production, and to predict if there are any possibilities of geo-hazard or reservoir failure. Calculation of the strength of hydrate bearing sediment is a complex process since the strength is a function of the hydrate saturation, strain rate, temperature, consolidation stress, grain size, density, and cage occupancy. The strength of hydrate bearing sediment is higher than the strength of sediment not containing hydrate and it is directly related to hydrate saturation [38]. At high hydrate saturations, the hydrate controls the strength and deformation characteristics in hydrate bearing sediments. Depressurization induced gas production from hydrate bearing sediments results in corresponding increase in effective stress and sediment compaction follow [39]. In this paper, we present our results for stress accumulations in reservoir through Figs. 7 and 8 in xx and yy directions respectively. These figures highlight total stress accumulations during 50 years of hydrate dissociation. Presence reduction, hydrate dissociation and fluid movement had direct impact on these stresses. The impact was more significant in the case of stresses in xx direction. Large accumulation of stress in hydrate layer in xx direction is clearly visible in Fig. 7. Accumulation of stress in yy direction

320

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

Fig. 4. Graphical representation of hydrate saturation (%) – after (a) 3 years (b) 10 years (c) 25 years and (d) 50 years.

stability, compaction or deformation of reservoir. Fig. 9 shows that changes in effective stress in xx direction (effective Sxx) took place in large areas of reservoir. Fig. 10 illustrates that changes in effective stresses in yy direction (effective Syy) were mostly confined in hydrate layer. Reduction in effective stresses in both directions indicates that compaction in hydrate layer happened during hydrate dissociation. Reduction in effective Sxx and effective Syy were around 4 MPa. Large increase in effective Syy near reservoir boundaries might be due to the presence of boundaries, which affects overall calculations. Strength of hydrate-filled sediments reported by Ebinuma et al. [41] suggests that this change in effective strength might be well within safe limits for any mechanical failure in the reservoir.

CH4 Production Rate (ST m3/day)

Time vs CH4 production rate

Time (years)

4.5. Comparison of results with other theoretical studies

Fig. 5. Time (years) vs gas production rate (ST m3/day).

was not much and it was mostly concentrated near reservoir boundaries as visible in Fig. 8. 4.4.1. Effective stresses To study geomechanics of the system, effective stress calculation has been implemented into RCB according to Terzaghi’s Principle [40]. According to this principle, effective stress controls the mechanical failure of rock and is defined as:

r0ij ¼ rij  P0 dij

ð14Þ

According to this definition, a tensile fracture will happen if the minimal effective stress is negative and its absolute value is greater than tensile strength of the formation [40]. RCB hence helps in making entire process faster without any delay in the study of reservoir failure criteria. Figs. 9 and 10 illustrate effective stresses in xx and yy directions respectively. Effective stress is needed in the studies of reservoir

While the main purpose of this paper have been to illustrate and discuss alternative implementation of thermodynamics and kinetic in non equilibrium hydrate filled sediments it could be useful to look at some other papers which report on similar studies using other simulators. For this purpose, we used the code comparison work of Anderson et al. [42] and compared our results with their results. Anderson et al. studied several possibilities of gas production from class 3 hydrate deposits for Mount Elbert, Alaska. One of their methods was purely depressurization based gas production from this low temperature hydrate reservoir. The study involved five different simulation codes, which were CMG STARS [43], HydrateResSim [44], MH-21 HYDRES [45], STOMP-HYD [46], and TOUGH+HYDRATE [47]. Some of the parameters for comparison purposes were volumetric gas production rates and lag times (time before significant gas production) for each code. Production rates reported by STOMP-HYD and HydrateResSim were less compare to CMG STARS, MH-21 HYDRES, and TOUGH+HYDRATE. Average lag time was 13.5 years for all these codes. Their study concluded that all codes showed remarkable agreement in the characteristics of the long term production simulations. RCB

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

321

Fig. 6. Graphical representation of gas pressure (MPa) – after (a) few hours (b) 1 month (c) 1 year (d) 10 years.

Fig. 7. Graphical representation of stress in xx direction Sxx (MPa) – (a) initially (b) after 50 years.

Fig. 8. Graphical representation of stress in yy direction Syy (MPa) – (a) initially (b) after 50 years.

produced promising results for both gas production rates as shown in Fig. 5 and lag time which was less than 3 years. Lag time is also depended on permeability models used in the code. Slight difference in permeability calculations may result in lag time reduction.

Similar theoretical studies of depressurization induced methane hydrate dissociation were also carried out by Moridis et al. [48] for Mount Elbert, Alaska. The hydrate deposits site was different than that used in the work of Anderson et al. [42]. TOUGH+HYDRATE

322

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

Fig. 9. Graphical representation of effective stress in xx direction Sxx (MPa) – (a) initially (b) after 50 years.

Fig. 10. Graphical representation of effective stress in yy direction Syy (MPa) – (a) initially (b) after 50 years.

[47] was used as hydrate simulator by Moridis et al. Conditions of reservoir site were unfavorable to hydrate dissociation and gas production as the reservoir temperature was very low (2.3– 2.6 °C). Main focus of the study was on the gas production rates. Evolution of temperature and pressure with time, change in hydrate and gas saturations in reservoir were also reported in this work. Final conclusions were that hydrate accumulations with low hydrate saturation appeared to be more desirable potential targets for a successful long term field test of production from colder, permafrost associated hydrates and initial deposit temperature was the most important factor determining production performance. Rutqvist et al. [49] studied class 3 hydrate deposits at Mallik and Mount Elbert deposits, Alaska. In their work, combination of TOUGH+HYDRATE simulator and FLAC3D [50] commercial geomechanical code were used. Simulation period was only 5 years which was relatively shorter. Low gas productions from both Mallik and Mount Elbert sites were reported. Their study also focused on geomechanical impacts of depressurization on reservoir and concluded that the depressurization of the hydrate reservoir results in vertical compaction of the reservoir and in increase of shear stress. The magnitude of vertical compaction and shear stress depend on the magnitude of depressurization and the elastic properties of the reservoir and overlying formations. One conclusion was that in both Mallik and Mount Elbert, the higher shear stress may lead to shear failure in the hydrate dissociation zone. Our results in Figs. 9 and 10 for effective stresses also confirm compaction in the hydrate layer during production. Results from RCB for stress accumulation shows that very high stress accumulation has occurred in hydrate layer during 50 years period, a trend which might lead to fracture or reservoir failure. A more detailed study of geomechanical stability of the reservoir, including measurements of tensile strengths is required in order to give qualified answers.

5. Conclusion Hydrates in porous media are not able to reach thermodynamic equilibrium and natural gas hydrate reservoirs are there generally in a stationary situation which depends on fluid flow, cap-rock of clay sealing and possible fracture patterns that support new hydrocarbon supply and/or fracture/fault systems which bring hydrate in contact with groundwater and corresponding local hydrate dissociation. We have proposed an alternative way of treating hydrates in reservoir simulations which makes use of well known logistics from reactive transport simulators. Kinetic rates for different hydrate dissociation ‘‘reactions’’ are derived from phase field theory simulations and similar for different hydrate formation ‘‘reactions’’. After comparing our results with the results of other hydrate simulators it has been demonstrated that our kinetic approach for hydrate description in RCB works very well. It is competent of handling hydrate kinetic ‘‘reactions’’ and is capable of generating all important results necessary for detailed analysis of a gas production scenario. Implicit geomechanical calculation makes RCB a good tool to predict reservoir stability during this gas production. In conclusion it has been demonstrated that after our extension of the RCB reservoir simulator is capable of functioning as a hydrate simulator on at least the same level of functionality as other academic hydrate simulators but with a more appropriate handling of thermodynamics and phase transition. At this stage it can handle simple hydrate reactions for hydrate formation and dissociation. Extensions to hydrate dissociation towards under saturated water and other relevant phase transitions are in progress and will successively be implemented. So far only kinetic rates for CH4 and CO2 hydrates are implemented but extensions are trivial since the implemented thermodynamic packages for fluid and hydrates are general and not limited to these two components. It does, however, require

A. Chejara et al. / Energy Conversion and Management 68 (2013) 313–323

development of kinetic rate equations from phase field theory simulations. Acknowledgement We acknowledge the grant and support from Research Council of Norway through the following Projects 178008/I30, 804831 and 801445. References [1] Makogon YF. Perspectives for the development of gas hydrate deposits. In: Fourth Canadian Permafrost Conference, Calgary, March 2–6, 1981. [2] Kvenvolden KA. A review of the geochemistry of methane in natural gas hydrate. Org Geochem 1995;23(11/12):97–1008. [3] Kvamme B, Kuznetsova T. Hydrate dissociation in chemical potential gradients: theory and simulations. Fluid Phase Equilibr 2004;217:217–26. [4] Graue A, Kvamme B, Baldwin B, Stevens J, Howard J, Aspenes E, et al. MRI visualization of spontaneous methane production from hydrates in sandstone core plugs when exposed to CO2. SPE J 2008;13(2):146–52. [5] Holder G, Angert P. Simulation of gas production from a reservoir containing both gas hydrates and free natural gas. Paper 11105-MS presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 26–29 September, 1982. [6] Burshears M, O’Brien TJ, Malone RD. A multi-phase, multi-dimensional, variable composition simulation of gas production from a conventional gas reservoir in contact with hydrates. Paper 15246-MS presented at the SPE Unconventional Gas Technology Symposium, Louisville, Kentucky, 18–21 May, 1986. [7] Yousif MH, Abass HH, Selim MS, Sloan ED. Experimental and theoretical investigation of methane-gas-hydrate dissociation in porous media. SPE Reservoir Eng 1991;6(1):69–76. [8] Kim HC, Bishnoi PR, Heidemann RA, Rizvi SSH. Kinetics of methane hydrate decomposition. Chem Eng Sci 1987;42(7):1645–53. [9] Moridis GJ. Numerical studies of gas production from methane hydrates. SPE J 2003;8(4):359–70. [10] Hong H, Pooladi-Darvish M, Bishnoi PR. Analytical modeling of gas production from hydrates in porous media. J Can Petrol Technol 2003;42(11):45–56. [11] Hong H, Pooladi-Darvish M. Simulation of depressurization for gas production from gas hydrate reservoirs. J Can Petrol Technol 2005;44(11). http:// dx.doi.org/10.2118/05-11-03. [12] Uddin M, Coombe D, Wright F. Modeling of CO2-hydrate formation in geological reservoirs by injection of CO2 gas. J Energ Resource Technol 2008;130(3):032502. [13] Methane Hydrate Reservoir Simulator Code Comparison Study. . [14] Saaltink M, Batlle M, Ayora C. RETRASO, A code for modeling reactive transport in saturated and unsaturated porous media. Geological Acta 2004;2(3):235–51. [15] Svandal A, Kuznetsova T, Kvamme B. Thermodynamic properties and phase transitions in the H2O/CO2/CH4 system. Fluid Phase Equilibr 2006;246(1– 2):177–84. [16] Kvamme B. Initiation and growth of hydrate from nucleation theory. Int J Offshore Polar Eng 2002;12(4):256–62. [17] Kvamme B. Droplets of dry ice and cold liquid CO2 for self transport to large depths. Int J Offshore Polar Eng 2003;13(1):1–8. [18] Kvamme B, Tanaka H. Thermodynamic stability of hydrates for ethylene, ethane and carbon dioxide. J Phys Chem 1995;99:7114–9. [19] Kvamme B, Kuznetsova T, Kivela PH. Adsorption of water and carbon dioxide on hematite and consequences for possible hydrate formation. Phys Chem Chem Phys 2012;14(13):4410–24. [20] Svandal A. PhD thesis submitted to University of Bergen, Bergen, Norway; 2006, p. 37. Available at Bergen Open Research Archieve (BORA-UiB) on . [21] Baig K. Master thesis submitted to University of Bergen, Bergen, Norway; 2009, p. 85. Available at Bergen Open Research Archieve (BORA-UiB) on . [22] M. Qasim, B. Kvamme, K. Baig, Modeling dissociation and reformation of methane and carbon dioxide hydrate using phase field theory with implicit hydrodynamics, Proceedings of the 7th International Conference on Gas Hydrates (ICGH 2011), Edinburgh, Scotland, United Kingdom, July 17–21, 2011.

323

[23] Qasim M, Kvamme B, Baig K. Phase field theory modeling of CH4/CO2 gas hydrates in gravity fields. Int J Geol 2011;5(2):48–52. [24] Svandal A, Kvamme B, Granasy L, Pusztai T, Buanes T, Hove J. The phase-field theory applied to CO2 and CH4 hydrate. J Cryst Growth 2006;287(2):486–90. [25] Tegze G, Pusztai T, Toth G, Granasy L, Svandal A, Buanes T, et al. Multiscale approach to CO2 hydrate formation in aqueous solution: Phase field theory and molecular dynamics. J Chem Phys 2006;124:234710. http://dx.doi.org/ 10.1063/1.2207138. [26] Saaltink MW, Benet I, Ayora C. RETRASO, Fortran code for solving 2D reactive transport of solutes. Users´s guide. Barcelona, E.T.S.I. Caminos, Canales y Puertos, Universitat Politecnica de Catalunya and Instituto de Ciencias de la Tierra, CSIC; 1997, p. 90. [27] Saaltink MW, Ayora C, Olivella S. RetrasoCodeBright (RCB), User’s guide, Barcelona, E.T.S.I. Caminos, Canales y Puertos, Universitat Politecnica de Catalunya and Instituto de Ciencias de la Tierra, CSIC; 2005. [28] Olivella S, Gens A, Carrera J. CodeBright User’s Guide. Barcelona, E.T.S.I. Caminos, Canales y Puertos, Universitat Politecnica de Catalunya and Instituto de Ciencias de la Tierra, CSIS; 1997. [29] Olivella S, Carrera J, Gens A, Alonso EE. Non-isothermal multiphase flow of brine and gas through saline media. Transport Porous Med 1994;15:271–93. [30] Olivella S, Gens A, Carrera J, Alonso EE. Numerical formulation for a simulator (CODE_BRIGHT) for the coupled análisis of salinemedia. Eng Comput 1996;13(7):87–112. [31] GiD user manual, CIMNE – International Center For Numerical Methods in Engineering, [email protected], Barcelona, Spain. ‘‘GiD – The personal preand post processor’’; 2012. . [32] Bear J. Dynamics of fluids in porous media. New York: American Elsevier; 1972. p. 132–5. [33] van Genuchten R. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 1980;44:892–8. [34] Horvath AL. Aqueous electrolyte solutions physical properties estimation and correlation methods, Ellis Holword, Chichester; 1985. [35] Hunter RB, Collett TS, Boswell RM, Anderson BJ, Digert SA, Pospisil G, et al. Mount elbert gas hydrate stratigraphic test well, Alaska North Slope: overview of scientific and technical program. J Mar Petrol Geol 2011;28(2):295–310. [36] Sloan ED. Clathrate hydrates of natural gases. NY: Marcel Dekker; 1998. [37] Alp D, Parlaktuna M, Moridis GJ. Gas production by depressurization from hypothetical Class 1G and Class 1W hydrate reservoirs. Energ Convers Manage 2007;48:1864–79. [38] Winters WJ, Pecher IA, Waite WF, Mason DH. Physical properties and rock physics models of sediment containing natural and laboratory-formed methane gas hydrate. Am Mineral 2004;89:1221–7. [39] Waite WF, Santamarina JC, Cortes DD, Dugan B, Espinoza DN, Germaine J, et al. Physical properties of hydrate-bearing sediments, Rev. Geophysics 47 (RG4003), p. 1–38. doi: 10.1029/2008RG000279. [40] Terzaghi K. Theoretical soil mechanics. New York: John Wiley and Son Inc.; 1943. p. 11–5. [41] Ebinuma T, Kamata Y, Minagawa H, Ohnuma R, Nagao J, Narita H. Mechanical properties of sandy sediment containing methane hydrate. In: Proc. 5th Int. Conf. on Gas Hydrate, Trondheim (Norway), vol. 3. 958–961 (Paper ref.3037), 13–16 June, 2005. [42] Anderson B, Kurihara M, White MD, Moridis GJ, Wilson SJ, Pooladi-Darvish M, et al. Regional long-term production modeling from a single well test, mount elbert gas hydrate stratigraphic test well, Alaska North Slope. Mar Petrol Geol 2011;28:493–501. [43] CMG – Computer Modeling Group: STARS User Manual, Calgary, Alberta, Canada; 2004. [44] Moridis GJ, Kowalsky MB, Pruess K. HydrateResSim Users Manual: A Numerical Simulator for Modeling the Behavior of Hydrates in Geologic Media. Department of Energy, Contract No. DEAC03-76SF00098. Lawrence Berkeley National Laboratory, Berkeley, CA; 2005. [45] MH-21 Research Consortium. . [46] White MD, Oostrom M. STOMP Subsurface Transport over Multiple Phases: User’s Guide PNNL-15782 (UC-2010). Pacific Northwest National Laboratory, Richland, Washington; 2006. [47] Moridis GJ, Kowalsky MB, Pruess K. TOUGH+HYDRATE v1.0 User’s Manual: A Code for the Simulation of System Behavior in Hydrate-Bearing Geologic Media. Report LBNL-00149E. Lawrence Berkeley National Laboratory, Berkeley, CA; 2008. [48] Moridis GJ, Silpngarmlert S, Reagan MT, Collett T, Zhang K. Gas production from a cold, stratigraphically-bounded gas hydrate deposit at the mount elbert gas hydrate stratigraphic test well, Alaska North Slope: implications of uncertainties. Mar Petrol Geol 2011;28:517–34. [49] Rutqvist J, Moridis GJ, Grover T, Collett T. Geomechanical response of permafrost-associated hydrate deposits to depressurization-induced gas production. J Petrol Sci Eng 2009;67:1–12. [50] FLAC3D, Itasca Consulting Group, Inc.. .