Simulations of Variable Bottomhole Pressure ... - ACS Publications

Mar 4, 2011 - Production from the Double-Unit Mount Elbert, Milne Point Unit, ... Kelly Rose,. † ... 2007 Mount Elbert field programs provide limite...
0 downloads 0 Views 4MB Size
ARTICLE pubs.acs.org/EF

Simulations of Variable Bottomhole Pressure Regimes to Improve Production from the Double-Unit Mount Elbert, Milne Point Unit, North Slope Alaska Hydrate Deposit Evgeniy M. Myshakin,*,†,‡ Brian J. Anderson,†,§ Kelly Rose,† and Ray Boswell† †

National Energy Technology Laboratory, 3610 Collins Ferry Road, Morgantown, West Virginia 26507, United States URS, 626 Cochrans Mill Road, Pittsburgh, Pennsylvania 15236, United States § Department of Chemical Engineering, West Virginia University, P.O. Box 6102, Morgantown, West Virginia 26506, United States ‡

bS Supporting Information ABSTRACT: Gas production was predicted from a reservoir model based on the Mount Elbert gas hydrate accumulation located on the Alaska North slope at various simulator submodels and production scenarios. Log, core, and fluid measurements were used to provide a comprehensive reservoir description. These data were incorporated with experimentally derived saturations, porosities, permeability values, parameters for capillary pressure, and relative permeability functions. The modeled reservoir exposed to depressurization at a constant bottomhole pressure (2.7 MPa) has shown limited production potential due to its low temperature profile. To improve production the bottomhole pressure was allowed to vary from 2.7 (above the quadruple point) to 2.0 MPa over a 15-year period. The results indicate that gas production was nearly doubled in comparison with a constant-pressure regime. Extensive ice formation and hydrate reformation that could severely hinder gas production were avoided in the variable-pressure regime system. A use of permeability variation coupled with porosity change is shown to be crucial to predict those phenomena at a reservoir scale.

1. INTRODUCTION Methane hydrate is a special class of gas inclusion compound formed with water molecules that is stable under high pressure and relatively low temperatures.1 It is abundant in nature in suboceanic sediments, primarily in continental shelf and slope regions, and in permafrost areas and is considered a potentially vast energy source.1,2 Methane hydrate formation and decomposition can be represented by the following reactions involving methane and water: CH4ðgasÞ þ nH2 OðliquidÞ T CH4 3 nH2 OðsolidÞ

ð1Þ

where n, the hydration number, is the number of water molecules per methane molecule (usually 6) in the hydrate lattice. Adaptations of three conventional natural gas production methods are usually considered as a means to recover the methane from hydrate-containing reservoirs.1 The first method is depressurization in which the pressure in the reservoir is lowered below the hydrate dissociation pressure. The second method is thermal stimulation in which the reservoir is heated above the dissociation temperature of the hydrate. The third method involves injection of inhibitors such as methanol or brine to shift the thermodynamic equilibrium and dissociate the hydrate. It is also thermodynamically possible to produce methane from hydrate by displacing it with a more favorable guest molecule such as carbon dioxide.3,4 This fourth, unconventional method also offers the potential for carbon sequestration.5 The depressurization production method has primarily been studied via long-term, reservoir-scale numerical simulations of r 2011 American Chemical Society

methane production from subsurface deposits using both constant-pressure and constant-mass-rate production methods.6-15 Two short-term field tests conducted during the 2002 Mallik and 2007 Mount Elbert field programs provide limited field data that are utilized by these models to better predict the effects of depressurization of natural gas hydrate-bearing sand reservoirs.11,13,16 Numerical simulations seek to emulate the conditions studied at field sites such as Mallik and Mount Elbert, but also theoretical relationships that vary the reservoir composition, architecture, and contents. As a result, these numerical simulations have shown that the constant-mass-rate method implying a withdrawal of fluids at a prescribed rate works well for hydrate formations with an underlying layer of free gas or water; however, without such a layer, the constant-pressure method can be used to achieve satisfactory production results. The advantage of the constant-pressure method is that the production rate can adjust itself to increased formation permeability that accompanies continuous hydrate phase decomposition.8,17 Also, maintaining a constant pressure at the well (Pw) above the quadruple point (PQ = 2.7 MPa) eliminates the possibility of ice formation during decomposition.8,17 Since the gas and water fluid rates out of the well are subject to timedependent phase mobilities and the pressure difference between the bottom of the well and the moving interface boundary, maximizing the Pw - Pinit (where Pinit is the initial reservoir Received: October 14, 2010 Revised: January 7, 2011 Published: March 04, 2011 1077

dx.doi.org/10.1021/ef101407b | Energy Fuels 2011, 25, 1077–1091

Energy & Fuels pressure) difference, ΔPw, usually provides effective depressurization of hydrate accumulations. On the other hand, field experiments conducted at the Mallik site and associated numerical simulations have shown that thermal stimulation by means of wellbore heating is an ineffective method for long-term gas production.9,18,19 This is primarily because of costly, wasted heating of the rock that contains the hydrate and diminishing effectiveness of the method as the dissociation interface propagates further from a wellbore. However, in recent reservoir simulations using depressurization, thermal stimulation has proven to be useful for reinstating gas production rates that have decreased or become zero because of hydrate reformation around the well.9,18,20 This hydrate reformation, also called secondary hydrate formation, was predicted in simulations utilizing a homogeneous reservoir characterization (uniform porosity, saturation, and permeability) and attributed to the endothermic nature of the hydrate dissociation reaction, compounded by additional cooling due to the Joule-Thomson effect accompanying the expanding gas as it nears and enters a well.9,18 In those simulations, the periodic addition of warm water increased production rate by dissociating the reformed hydrate barrier. There are several known simulators that have been developed and used for modeling gas hydrate production. Researchers at The Pacific Northwest National Laboratory and the Petroleum Engineering Department at the University of Alaska Fairbanks have adopted the multiphase simulator STOMP to extend its capabilities to consider gas hydrates (STOMP-HYD).21 The MH-21 Hydrate Reservoir Simulator (MH-21 HYDRES)22 has been developed by the National Institute of Advanced Industrial Science and Technology, Japan Oil Engineering Co., Ltd. and the University of Tokyo to study production from gas hydrate deposits. Another group studying Alaska North Slope gas hydrate resource potential as a part of the BP Exploration Alaska, Inc. (BPXA) research project in collaboration with the U.S. Department of Energy (DOE) has extended work begun at the University of Calgary and the University of Alaska Fairbanks to apply a commercially available simulator, CMG STARS,23 to model production from gas hydrate-bearing reservoirs.24 National Energy Technology Laboratory (NETL) has released a freeware, open-source, version of a code under the name HydrateResSim, which is available at the NETL Web site.26 This code is based on an earlier version of the TOUGHþHydrate code25 developed at Lawrence Berkeley National Laboratory (LBNL). There are differences in the capabilities and applications of these models, and over the past five years NETL has led an ongoing international effort to compare the predictions of these simulators.27,28 The numerical studies in this paper were carried out using a recent parallel version of the TOUGHþHYDRATE code25 to model the nonisothermal gas release, phase behavior, and flow of fluids and heat in complex geological media. The TOUGHþHYDRATE code can be used to simulate production from both natural gas hydrate deposits (reservoir simulations) and laboratory experiments of gas hydrate dissociation/formation processes in porous media (reactor simulations). The simulator accounts for heat and up to four mass components (water, methane, hydrate, and water-soluble inhibitors) divided among four possible phases (gas, aqueous, ice, and hydrate). The corresponding phase diagram of the water-methane-hydrate system is given in Figure S1 (Supporting Information). The model can describe a combination of three gas hydrate dissociation

ARTICLE

mechanisms (depressurization, thermal stimulation, and inhibitorinduced effects). The code involves a fully implicit formulation, continuous property updating, and Newton-Raphson iteration for the solution of the nonlinear coupled algebraic equations of fluid and heat flow. An interested reader can find concepts, underlying physics, governing equations, and details of the code implementation in the comprehensive manual available online.25 Although numerical simulations provide a low-cost way of studying gas production from hydrate-bearing formations, validation of the simulation results is often problematic considering the current lack of longer-term production data from hydratecontaining reservoirs9 and scarce field data. As noted above, only four short-term field-scale tests have been conducted in recent years, including tests at the Mallik research site in Northwest Canada in 2002, 2007, and 200811,29-31 as well as pressure transient tests at the Mount Elbert site in northern Alaska in 2007.32 As more field data characterizing actual gas hydrate deposits become available, more sophisticated reservoir models can be developed. In particular, availability of log porosity and hydrate saturation data provides a way to create heterogeneous reservoir models accounting for actual variability in porosity and hydrate saturation. In the literature there are a few works explicitly incorporating heterogeneity in the reservoir model. Recently, Anderson et al.13,15 compared production forecasts of gas hydrate reservoirs from the Alaska North Slope using heterogeneous porosity, saturation, and permeability descriptions with those created using homogeneous descriptions. These studies suggested greatly accelerated and enhanced gas production as compared to the prior studies. Kurihara et al.10,12 included heterogeneity with a different scale such as a fault, dip, and change in permeability for prediction of gas production from the Eastern Nankai Trough methane hydrate reservoirs. It was found that sealing capability, conductivity of faults, and formation dip are the factors significantly affecting short-term production. Reagan et al.17,33 investigated the impact of heterogeneity (both vertical and lateral) and found that vertical hydrate saturation variability derived from well log data results in enhanced gas production forecasts due to effects of small zones of greater in situ permeability. The simulations using random initialization or a geostatistical model show complex physical behavior and display poorer performance in comparison with a base homogeneous case. In this study, we model gas hydrate production performance for accumulations based on two subpermafrost hydrate deposits in the Alaska North Slope’s Milne Point Unit. In 2007, the U.S. DOE, BP Exploration (Alaska), and the U.S. Geological Survey collected open-hole formation pressure response data, well logs, cores, and pore water samples from the “Mount Elbert” stratigraphic test well.34 Analysis of the wireline log data confirmed that two hydrate-bearing sand reservoirs (units D and C) exist within the Sagavanirktok formation in that well. An extensive data set including hydrate saturation, water content, porosity, density, and permeability obtained from the well logs and cores34-36 was used to create a detailed geologic model for the numerical simulations of the reservoir. The Mount Elbert program also produced detailed characterizations of the litholostratigraphy,37 pore water geochemistry,38 gas geochemistry and composition,39 and the sediment physical properties.40 Limited production potential of the unit D or C hydrate zones were predicted from initial constant-P depressurization simulations, with the relatively low temperature of the reservoir and surrounding strata (2.3-3.9 °C) interpreted to be the primary 1078

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

ARTICLE

Table 1. Reservoir Properties, Initial Conditions, and Pertinent Model Parameters parameter

value

hydrate zone thickness combined

26.11 m

gas composition

100% CH4

water salinity (mass fraction)

0.005

initial hydrate saturation

heterogeneous derived from log data

intrinsic permeability of unit D and C, k

1100 mD horizontal (600 mD vertical) and based on perm-porosity correlation

porosity of unit D and C, φ

heterogeneous derived from log data

pore compressibility, Rp

6  10-9 Pa-1

intrinsic permeability of the intermediate shale layer, k porosity of the intermediate shale layer, φ

20 mD horizontal, 10 mD vertical 0.29 (averaged)

dry thermal conductivity, kdry

0.5 W/m K

wet thermal conductivity, kwet

3.1 W/m K

composite thermal conductivity,42 kΘ

pffiffiffiffiffi pffiffiffiffiffi kΘ ¼ kdry þ ð SA þ SH Þðkwet - kdry Þ þ φSI λI Pcap ¼ - P½ðSÞ - 1=λ - 11 - λ

capillary pressure model43

ðSA - SirrA Þ S ¼ ðSmaxA - SirrA Þ SirrA

heterogeneous derived from log data

λ

0.7743713

SmaxA

1

Pmax

104 Pa

1/P

1.1  10-3 Pa



modified Stone 3-phase model44



SA ¼

ðSA - SirrA Þ  ðSG - SirrG Þ ;S ¼ ð1 - SirrA Þ G ð1 - SirrG Þ

n

nA = 4.5013; nG = 3.1613

SirrG

0.02

SirrA

heterogeneous derived from log data

constraint on the gas production rate.27,41 The goal of this study is to investigate gas production from the hydrate reservoir using detailed geologic descriptions under production scenarios designed to enhance reservoir performance. These simulations utilize the close proximity of the hydrate units D and C and evaluate the effect of decreasing bottomhole pressure below the quadruple point on key characteristics of the reservoir during production. This paper is organized as follows. In the next section we describe details of the mesh discretization, initial conditions, and parameters. Then, an approach to derive porosity-permeability relationship is considered. The Hydraulic and Wettability Models subsection provides details on how effect of solid phase change in pore space is treated in respect to relative permeability and capillary pressure. Section 3 discusses results of gas production from the reservoir model at P/T conditions existing at the Mt. Elbert hydrate deposit site. The model utilizes either the derived porosity-permeability dependence or uniform permeability values. This is followed by results of the simulations performed at P/T conditions relevant to the L-Pad hydrate accumulation site at Prudhoe Bay Unit (PBU) on the North Slope Alaska south of the Mount Elbert site. The section ends with economic assessment of productivity performance at both constant-P and variable-P depressurization regimes. A summary of our findings and implications of the variable regime to gas production from hydrate reservoirs is given in Section 4.



krA ¼ ðSA Þn ; krG ¼ ðSG Þn

relative permeability

2. SIMULATION DETAILS 2.1. Reference Geologic System. The reference geologic system built for this study was based on log, core, and seismicbased data from the BP-DOE-USGS “Mount Elbert” Stratigraphic Test Well.34-36,40 The Mount Elbert well included the acquisition of 131.7 m (432 feet) of core and an extensive suite of wireline log data. Analysis of the core data and log-based interpretations confirmed that 21.3 m (70 feet) of thinly interbedded silt-rich sediments separates the upper 14.0 m (46 feet) thick hydratebearing unit D sand-rich reservoir from the 16.2 m (53 feet) thick sand-rich unit C. Gas hydrate was dominantly observed in the more sand-rich intervals with variable pore saturations up to 80%.35 There was no free gas evident within the system. 2.2. Description of Gas Hydrate Reservoir Model. Information on reservoir properties, initial conditions, and pertinent model parameters are given in Table 1. The parameters for relative permeability and capillary pressure functions were determined through history matching of a mutli-stage well test from the Mount Elbert stratigraphic test well using Schlumberger’s Modular Dynamics Formation Tester (MDT) wire line tool.13 Figure S2a and S2b in the Supporting Information shows relative permeabilities vs water saturation and capillary pressure vs water saturation dependencies, respectively. The gas hydrate deposit modeled here is based on the log and core data from the Mount Elbert site, and closely mimics those deposits in terms of depth, thickness, and petrophysical properties. For these simulations, 1079

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

Figure 1. Schematic of Mount Elbert hydrate deposit.

however, the boundary conditions are simplified to include impermeable shale units both above and below the accumulations. The reservoirs are modeled as a cylindrical accumulation including 14.07 m of unit D, 12.04 m of unit C, and 21.25 m of a semipermeable silt layer placed between the units (Figure 1). The inner radius of the system was the well radius equal to 0.1 m, and the outer radius was set at 550 m providing no flow and no heat exchange through vertical sides of the domain. The model has used a 2D representation taking advantage of the cylindrical symmetry, so the discretization of the domain into 152  239 = 36328 (r,z) grid blocks has resulted in 72 265 connections and 181 640 coupled nonlinear algebraic equations solved simultaneously at every Newton-Raphson iteration by the Biconjugate Gradient Squared solver, which provides robust convergence for the simulations performed in this work. Detailed information about implemented solvers and their advantages and disadvantages can be found in the manual25 and other publications.45,46 Units D and C are discretized into a 60-layer hydrate block with a 0.2345-m-thick individual layer and a 50-layer hydrate block with a 0.2404-m-thick individual layer, respectively. The intermediate silt-rich layer between the units includes an 85-layer block with a 0.25-m-thick individual layer. The radial discretization contains 152 grid blocks with logarithmically distributed lengths, so very fine meshing is provided along a wellbore. The model includes overburden and underburden blocks allowing heat, but no mass exchange with the hydrate units. The vertical discretization of the overburden and the underburden consists of 22 layers of different thickness (1  5 m; 6  2.5 m; 2  2.0 m; 3  1.0 m; 3  0.5 m; 6  0.25 m) with the last six 0.25-m-thick layers closest to the hydrate layers. So, both bounding units have 30-m thickness, which earlier calculations proved to be sufficient to provide heat exchange with a hydrate-bearing layer during a 30-year production period.7 The wellbore located at the center of the cylindrical domain was constructed vertically through unit D, the intermediate silt block, and unit C (Figure 1). The wellbore was represented using the pseudoporous-medium approach18 providing deviations in results no more than 5% from the

ARTICLE

accurate solution of the Navier-Stokes equation. Employing the equilibrium model25 and the depressurization method to induce hydrate dissociation, the topmost grid block (which is located just above the boundary line between the overburden and the upper unit D block) of the subdomain representing the wellbore was set as an internal boundary maintained at constant or variable bottomhole pressure (Pw). 2.3. Initial Conditions. The pressure in the sediment subsurface was assumed to follow a hydrostatic pore pressure distribution. This assumption is supported by measurements taken in natural hydrate deposits.47 To determine thermal distribution throughout the reservoir model the local geothermal gradient equal to 0.029 °C/m was used.34,48 For unit D this resulted in average pressure and temperature of 6.87 MPa and 2.79 °C, respectively, and 7.21 MPa and 3.81 °C for unit C, respectively. For the layers of the reservoir model proper initialization of initial P/T conditions (pressure and temperature gradients throughout the vertical dimension of the domain) is needed to achieve hydraulic, thermal, thermodynamic, and chemical equilibrium and ensure correct location of the layers relative to the base of gas hydrate stability. To do that a procedure similar to that reported in ref 7 was used. Pore water salinity was assumed to be 5 ppt based on the measurements made at the Mount Elbert well.38 The temperature depression induced by salt presence in the aqueous phase on the equilibrium pressure-temperature (P/T) relationship was computed internally in the code through the equation of Dickens and Quinby-Hunt.49 The model uses vertical heterogeneity of porosity, hydrate saturation (SH), and irreducible water saturation (SirrA) derived from well log data. This provides “layered” heterogeneity with the same parameters extended horizontally for the mesh layers inside the hydrate-bearing units D and C. This is still an approximation to actual 3D spatial variability of porosity/saturation of hydrate-bearing media, however, it is a more robust one than a homogeneous approach implying uniform properties of a hydrate accumulation. To the best of our knowledge the reservoir model used in this study is the most elaborate one utilizing well log data and both hydrate units. Recently Regan et al.33 used heterogeneous hydrate saturation of unit D to assess an effect of heterogeneity on gas production and Anderson et al.14,15 employed vertical variability of saturation, porosity, and irreducible water saturation for unit C hydrate deposit within the code comparison project activity.27,28 A direct use of well log data in the reservoir model would result in a very dense mesh, which makes simulations computationally prohibited. This dictates a necessity to average data to assign heterogeneous porosity, SH, and SirrA values to the individual layers of the hydrate-bearing units D and C. To that end well log data were interpolated with analytical functions representing polynomials (r2 > 0.99). Integration of the polynomial functions defined on intervals corresponding to a specific layer within units D and C provides a set of averaged values for the layers along the vertical dimension (Figure 2). Then the model was extended horizontally along the layers, such that those values were constant along that direction. For every layer the model data together with permeability values (see Section 2.4 below) are collected in a table and can be found in the Supporting Information. The resulted initial hydrate saturation and porosity distribution of the reservoir model are shown in Figure 3. 2.4. Permeability-Porosity Model. One of the key aspects affecting gas production from a reservoir is a relationship 1080

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

ARTICLE

Figure 3. Initial hydrate saturation (a) and porosity distribution (b) used in the reservoir model (over- and underburden are not shown).

Figure 2. Log and modeled porosities (φ), hydrate saturations (SH), irreducible (SirrA), and free water saturations across the vertical dimension of units D (a) and C (b).

between intrinsic permeability (permeability in the absence of gas hydrate) and porosity. Since our reservoir model incorporates heterogeneous porosity it is natural to relate different pore volumes available for flow with a change in permeability. To do so permeability values obtained with a core probe minipermeameter40 were correlated with logderived porosities measured at the same depth (as could best be determined). A correlation between core-measured and well log derived porosities is strong with porosity values determined from core plugs typically similar or slightly higher than the log values.40 It was found that in order to improve the permeability-porosity correlation, only porosity values between 0.3 and 0.42 were included in the analysis. Since 90% of the total log-derived porosity values for units D and C fall within that range this assumption seems justified. This porositypermeability correlation is depicted in Figure 4 together with an equation resulted from a linear regression analysis of the data set. The obtained equation was further employed to introduce vertical variability of intrinsic permeability for the hydrate units using corresponding model porosity values (Figure 2a and 2b).

According to the deduced equation (Figure 4), a 4% volume increase to a porosity value changed permeability by 60% on average. This shows rather a mild dependence of permeability on porosity change.1 Winters et al.40 reported a physical property analysis of samples from core section intervals covering the entire log depth. Permeabilities of the samples were determined by direct measurement (at net confining stress (NCS)) and by calculations using the method based on laser-grain-size analysis, LGSA.50 It was concluded that a 4% change in porosity altered permeability by an order of magnitude on average. A possible reason of deviation in the permeability-porosity dependency used in this work and reported by Winters et al. can be related to a selection of samples for the analysis. Figure 5 depicts histograms of sample porosities used in the NCS and LGSA analysis,40 those employed in a correlation with miniperm measurements (this work) and log porosities for units D and C.34,35 Porosities measured for the NCS/LGSA analysis fall at the tails of a normal distribution of the log porosity population and outside of one standard deviation of its mean. At those porosity values, evaluated permeabilities experience the most drastic changes resulting in a large sensitivity of permeability on porosity change. On the contrary, the deduced correlation using minipermeameter data was based on a larger set of porosity values that are within a standard deviation of the log porosity mean (Figure 5). 2.5. Hydraulic and Wettability Models. Lowering bottomhole pressure below the quadruple point could cause ice formation during hydrate decomposition because according to the phase diagram of the water-methane-hydrate system51 the equilibrium P/T conditions exist below the freezing point. Ice formation in hydrate-bearing porous media could lead to an 1081

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

Figure 4. Permeability measured with a core probe minipermeameter vs log porosity.

ARTICLE

of the pore space is distributed between gas (aqueous) and solid phases. Permeability reduction for the mobile phases occurs due to the presence of solid phases, (SA þ SG þ SH þ SI = 1), where SH and SI are hydrate and ice saturation, respectively. A shortcoming of the OPM model is that solid cannot be simultaneously present only in liquid or gas-filled pore space because according to the model liquid phase flow behaves as though solids deposition occurs entirely in what would otherwise be gas-filled pore space, while gas phase flow behaves as though solids deposition occurs entirely in what would otherwise be liquid-filled pore space. However, in spite of the limitation of the model, it is a popular choice because of its simplicity, which is useful when field data are insufficient. In the Evolving Porous Medium model (EPM) intrinsic permeability is modified in response to hydrate saturation change since changing saturation of solid phases alters porosity and permeability with time. The relative permeability function remains the same as for the OPM model except that aqueous (gas) saturation is adjusted according to: SAðGÞ ¼

SAðGÞ SA þ SG

ð2Þ

which reflects that now (effective) pore space is divided between liquid and gas phases only (thus, this model is free from the OPM model limitation). The relationship describing a change in permeability in response to a porosity change is based on solid precipitation studies:53,54   k φ - φc n ¼ ð3Þ k0 φ0 - φc

Figure 5. Histograms of porosity at Log (unit D and C), NCS/LSDA,40 and minipermeameter (this work) data sets.

adverse effect on gas production because of decreased permeability. On the other hand, exposure of hydrate-bearing sediments to lower pressure might result in higher gas production rates owing to a larger pressure drop between a wellbore and a dissociating interface. Another process that leads to solid phase emerging during gas production is known as secondary hydrate formation18,9 caused predominantly by the endothermic hydrate decomposition reaction and the Joule-Thomson effect of exiting gas. Salinity drop due to dilution by fresh water from dissociating hydrate can also contribute to the occurrence of hydrate reformation.8 Because of a strong impact of solid phase evolution in pore space on gas production it is important to consider how the simulator algorithms controlling hydraulic behavior and capillary pressure treat solid saturation change. In this work, the modified Stone 3-phase function was used to calculate relative permeabilities with parameters given in Table 1. In the simulations this function was used in connection with two models. The Original Porous Medium model (OPM) treats solid phases as other mobile phases except that solids are stationary. As a result, solid phases are accounted as a relative permeability correction in relation to effective permeability of mobile phases.25 The relative permeability depends only on aqueous, SA, (gas, SG) saturations and is independent of how the remaining fraction (1 - SA (SG))

where k0 and φ0 are intrinsic permeability and porosity of hydrate-free, and k and φ are intrinsic permeability and porosity of hydrate filled porous media; φc is the critical porosity at which permeability becomes zero. A 0.05 value used for φc was derived from the concept of pore bodies and pore throats.55 As gas hydrate continues to precipitate in the interstitial space it completely clogs the pore throats and disconnects the porous network while there is still space available in the pore bodies. The n value can vary significantly from 2 to 3 (mild dependence)52 to 10 or more56 and is an empirically derived parameter. To estimate n, the unit D and C density porosities (φ0) and NMR porosities (φ) from the downhole logging measurements were used.34,35 Those data provide a robust technique to evaluate hydrate content in pore space. The related permeability data represent the in situ NMR readings of KSDR (SchlumbergerDoll Research) permeability (k) and permeability obtained with a core probe minipermeameter (k0) taken at the same depth (accurate core to log depth control could be an issue here). It is evidently that differences between intrinsic permeability values measured by a minipermeameter and in situ (in the presence of gas hydrate) KSDR permeability are primarily the result of gas hydrate presence. By substituting those porosity and permeability data into eq 3 the n parameter was determined to be 5.67 ( 1.86, which denotes a moderate intrinsic permeability adjustment upon hydrate deposition in porous media. The capillary pressure model used in the simulations employs the van Genuchten function with parameters deduced from experiment13,15 (Table 1). The presence of solid (hydrate þ ice) phases in the pore space changes capillary suction and requires corrections to capillary pressure in hydrate and ice free porous media. In the simulations the capillary pressure scaling 1082

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

ARTICLE

was done by the Leverett method57 using the following relationship: sffiffiffiffiffiffiffi k0 φ Pcap, 0 ð4Þ Pcap ðSA Þ ¼ kφ0 where Pcap,0 is the capillary pressure computed with the van Genuchten function for ice and hydrate free medium and corresponding to k0 and φ0. The k0/k ratio is calculated using eq 3 assuming that the effect of ice formation on intrinsic permeability mimics that of hydrate. This assumption is based on the results of experiments reporting that hydrate formation is primarily pore-filling, not grain-cementing58-62 and ice forming convex shapes within pore centers that minimize surface area.58,63 2.6. Production Cases. To induce methane hydrate decomposition, the bottomhole pressure (Pw) was set at 2.7 MPa, just above the quadruple point to eliminate possible ice formation and maximize pressure drop at a wellbore. At that constant pressure value the simulations were run over a 30-year period corresponding to usual life-span of a wellbore and represent the “0” scenario. To explore a reservoir response to further pressure decline at a wellbore, two new scenarios were designed. The first one uses the same constant pressure regimes as for the “0” scenario for the first 5 years, after which bottomhole pressure linearly decreases from 2.7 to 2.0 MPa over the next 15 years, then stays constant to the end of a simulation (30 years). This scenario is denoted as “5”. The second one denoted “10” is similar to “5” except that a variable pressure is employed later at a 10-year simulation time point and lasts until 25 years of simulation time. After that the pressure stays constant again until the end of a simulation. Further in the text the “5” and “10” scenarios are referred to as scenarios with variable bottomhole pressure regimes. Depending on how effective permeability of mobile phases and capillary pressure are treated three methods were employed in the simulations. The first one (“OR”) uses the OPM model without capillary pressure scaling. The second one (“EV”) uses the EPM model without capillary pressure scaling. The final method (“EVS”) uses the EPM model with capillary pressure scaling included. Detailed descriptions of the models are given in Section 2.5. Combinations of these methods with the designated scenarios provide nine production cases: OR0, OR5, OR10; EV0, EV5, EV10; EVS0, EVS5, EVS10.

3. RESULTS AND DISCUSSION Figures 6 and 7 show production rates and cumulative volumes of gas produced at a wellbore for the cases of interest. It is immediately evident that lowering bottomhole pressure values below the quadruple point turned out to be beneficial for the production at all cases using the “5” and “10” scenarios. After both 5-year and 10-year time points a linear decrease of the pressure from 2.7 to 2.0 MPa leads to substantial increase in production rate and total volume of produced gas. The “5” scenario shows greater gas rate and total volume than the “10” one. Decreasing bottomhole pressure below the quadruple point also resulted in a higher gas/water ratio. The cumulative volumes of gas produced over a 30-year period are given in Table 2. In comparison with the OR0/EV0/EVS0 cases the gas volume has increased 1.59/1.60/1.53 and 1.73/1.88/1.76 times at OR10/ EV10/EVS10 and OR5/EV5/EVS5 cases, respectively (Table 2).

Figure 6. Rates of gas produced at a wellbore at (a) OR0(5)(10), (b) EV0(5)(10), and (c) EVS0(5)(10) cases.

The reported volume gains are substantial taking into account that for the “5” and “10” scenarios the hydrate-bearing deposit was exposed to slowly declining bottomhole pressure stretched over a 15-year period. In absolute values the reported cumulative volumes of gas noticeably deviate for different methods (OR, EV, and EVS) engaging the same scenarios; “0”, “5”, and “10” (Table 2). The OR method using the original porous medium model, which is a “default” choice for reservoir simulations in the absence of field data, predicts the highest rates. The EV method using the evolving porous medium model without capillary pressure scaling demonstrates the lowest rate values. If scaling is applied to accommodate a change in porosity (the EVS method) the rates fall between the EV and OR results. The OR method provides an upper estimation of production potential over EV and EVS because relative permeability correction due to hydrate presence at the OPM model leads to higher effective permeability of mobile phases than it is computed with absolute permeability adjustment at the EPM one. For all cases, an initial production period is characterized by low gas production rates in approximately the first 5 years (Figure 6). During that period the flow out of a wellbore is dominated by liquid water fraction. The gas rate starts to pick up 1083

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

ARTICLE

Figure 7. Cumulative volume of gas produced at a wellbore at (a) OR0(5)(10), (b) EV0(5)(10), and (c) EVS0(5)(10) cases.

Table 2. Cumulative Volumes of Gas Produced during a 30Year Period for Studied Cases OR0 OR10 OR5 8

EV0 EV10 EV5 EVS0 EVS10 EVS5

3

10 m 0.7989 1.2693 1.3800 0.3802 0.6060 0.7132 0.5680 0.8664 1.0014

after pronounced horizontal interfaces are developed at the top and the bottom of the reservoir in response to heat flow from surrounding strata. Appearance of such interfaces is a distinctive feature of Class 3 reservoirs characterized by a single hydratebearing layer sandwiched between shale layers and without mobile gas or water interval underneath.8,9 In Figure 8, pressure, hydrate, and gas saturation distributions are given for the EVS0/ 5/10 cases after 15 years of production. To trace evolution of those saturations and pressure, Figure S3 (Supporting Information) depicts their distributions in the reservoir after 20 years of production. Hydrate saturation distribution revealed that in unit D hydrate dissociation resulted in a horizontal interface distinctly evolved at the top of the reservoir while in the bottom of unit C hydrate remained less decomposed because of smaller pressure drop inside the unit (Figure 8, third row). Moreover, the hydrate saturation distributions show preferential

dissociation of select layers with high permeability that provides increased surface area for delivery of pressure perturbations within the reservoir. This hydrate dissociation pattern is attributed to the heterogeneous porosity and hydrate saturation where layers possessing low saturation and high permeability values reach complete hydrate dissociation faster than in the adjacent areas.15 The middle row of the Figure 8 displays gas saturation distributions over the first 60 m from a wellbore for the EVS0, EVS10, and EVS5 cases. It is noticeable that methane is accumulated within select layers. However, gas saturation level is relatively low inside the reservoir and gas released due to the dissociation reaction is readily delivered to a wellbore for production through these layers. The average recovery values computed as a ratio between volumes of gas produced at a wellbore and released inside the reactor due to hydrate decomposition are consistently within a high 0.8-0.9 range during the production period. Because of high recovery there is relatively low gas accumulation inside the reservoir during a 30-year production period. The only area of relatively high (0.5) gas concentration spreading over several meters is developed for the EVS5 case (Figure 8) showing the highest production rate among those three cases. Over the next 15 years of production that high gas saturation area extends only slightly. The low level of gas accumulation inside the reservoir found in this study is in contrast to predictions made for early Class 3 reservoir simulations employing homogeneous representations of porosity/ hydrate saturation and displayed higher production rates.6,8,17,64 In those studies high gas saturation is usually observed at the top of a reservoir due to the buoyancy forces and consequently becomes a main source for gas produced at a wellbore. High rates of gas flow and that associated with the Joule-Thomson cooling promote initial secondary hydrate formation near a wellbore. Yet such high gas accumulation at the top of a reservoir is possible due to impermeable boundary conditions between the hydratebearing media and the overburden. This is an approximation because in nature gas would inevitably migrate into surrounding strata and reduce the production rate. Its further fate is unknown and might be a subject of a separate investigation. In this regard the predictions made in this paper are less influenced by the impermeable boundary condition approximation because of low level of gas accumulation at the top of the reservoir. The third row of Figure 8 captures pressure distribution at the EVS0, EVS10, and EVS5 cases. At EVS0, the depressurization has propagated effectively in the upper unit D where a pressure gradient only slightly increases above 2.7 MPa as traced away from a wellbore. In unit C, pressure remains within 1 MPa above the bottomhole pressure after 15 years of production. The reason for that is a presence of a 21.25-m low permeability silt unit which acts as a hydraulic barrier for effective depressurization of the lower hydrate unit. In a typical Class 3 reservoir sandwiched between over- and underburden the largest pressure drop is induced at the top of a well near the topmost grid block. After beginning of production a dynamic pressure distribution is established throughout the rest of the grid blocks representing a wellbore completion. As dissociation continues the permeability of media around a well changes from very low values (due to hydrate presence) to high ones. This process provides more effective depressurization of dissociating hydrate with an interface receding from a wellbore. Additionally, the development of the horizontal interfaces8,9 (due to heat flow from over- and underburden) create highly permeable, free of hydrate areas assisting fast depressurization of far hydrate regions of a reservoir. 1084

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

ARTICLE

Figure 8. Hydrate saturation, gas saturation, and pressure distribution (Pa) for the EVS model after 15 years of production at three different scenarios.

In the model the intermediate silt layer always stays at very low permeability and, therefore, it hinders effective depressurization via porous media. So, unit C is exposed to smaller pressure drop as provided by the lower portion of a wellbore (where initial pressure of a well grid block adjacent to the top of unit C is at 3.6 MPa). This resulted in slower decomposition of unit C than in unit D (Figure 8, the upper row). Slow hydrate dissociation is also aggravated by lower intrinsic permeability of unit C. Among the two units the C one possesses almost 50% smaller average permeability (454 vs 819 mD), which reflects smaller on average porosity value (Figure 2a and 2b). It is also worth noting that throughout almost the entire thickness of unit C, free water saturation oscillates within a narrow [0.05-0.1] range (Figure 2b), while for unit D those values fall into a wider [0.1-0.3] bracket (Figure 2a). As a result low free water level hinders depressurization in unit C. In the EVS10 and EVS5 cases, pressure distributions given at the same scale as EVS0 are presented in Figure 8. The pressure change is expectably the lowest for the EVS5 case, for which the bottomhole pressure slowly decreased starting from a 5-year mark (attempts to lower pressure earlier do not improve production because of the limited area of dissociating interface and ice formation emergence causing a decrease in a gas rate). Noticeably, throughout unit C pressure still remains slightly above the quadruple point and stays above that level during the whole 30-year simulation period. Thus, only the upper unit D was exposed to lower than 2.7 MPa pressure value at which ice formation is possible. A thorough analysis of simulation data

reveals absence of extensive ice formation in unit D during gas production at the OR5/10, EV5/10, and EVS5/10 cases. This confirms that slow decrease of the bottomhole pressure from 2.7 to 2.0 MPa was beneficial for production performance at the “5” and “10” scenarios. We further investigated this effect by imposing a more steep pressure change (from 2.7 to 1.0 MPa) over the same time period after 5- and 10-year marks in an attempt to increase gas production while keeping the rest of the parameters and initial conditions the same. The results have shown that after several years of production at the bottomhole pressure reaching 1.4 MPa, ice and secondary hydrate formation caused a halt of production. Thus, the variable bottomhole pressure regime used at the “5” and “10” scenarios seems to be an optimal choice to avoid ice and secondary hydrate formation and benefit from higher production rates. 3.1. Simulations with Enhanced Permeability. Permeability of a hydrate-bearing reservoir is a key parameter controlling production performance.65,66 The heterogeneous permeability used in the reservoir model provided average values equal to 819 and 454 mD for units D and C, respectively. Those values represent a rather conservative estimate of average permeability. For example, Winters et al.40 reported average “plug” intrinsic permeability to nitrogen values of 1700 mD and 675 mD for units D and C, respectively. To explore production response to increased permeability, uniform intrinsic permeability values, equal to 1100 mD in the horizontal direction and 600 mD in the vertical one, were assigned to units D and C. Those numbers are typical evaluations of intrinsic permeability for coarse-grained 1085

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

ARTICLE

Figure 9. Hydrate (a), gas (b), and ice (c) saturations and temperature distribution (d) in unit D after 15 years of production at the OR10i case.

sand and often taken as “default” idealized values when permeability measurements are not available.7,8,13 The original porous medium model employed in the OR0/5/10 cases was used in these simulations because OPM is a common choice when field data allowing determination of permeability dependence upon hydrate presence in pore space are absent. The new cases are designated as OR0i, OR5i, and OR10i. A comparison of production rates computed at the OR0/5/10 and OR0i/5i/10i cases has shown almost a 2-fold increase in production rates at each “0”, “10”, and “5” scenario. The results also reveal wide formation of ice in unit D after about 11 and 15 years of production at the OR5i and OR10i cases accordingly. Figure 9 depicts hydrate (a), gas (b), and ice (c) saturation distribution, and temperature distribution (d) along the first 60 m of unit D after 15 years of production at the OR10i case. The heterogeneous porosity and hydrate saturation promote a “layered” pattern of hydrate dissociation that creates layers filled with gas interspersed within hydrate-bearing media (Figure 9a and 9b). In the middle of unit D, a long multimeter ice block is formed in the region free of hydrate (Figure 9c) in response to temperature falls below 0 °C in that area (Figure 9d). An analysis of the results showed that ice first formed after 15 years of production and further stretched over tens of meters following the receding dissociation interface. In spite of significant extent, ice formed in unit D does not block the top of the reservoir where gas predominantly flows to a wellbore (Figure 9c, the arrow). Consequently, at the OR10i case the ice formation pattern provides a possibility for mobile phases to be produced at a wellbore with little impediment. The similar ice evolution was

predicted at the OR5i case except that ice started to form several years earlier and additionally developed around a wellbore in unit D after 15 years of production. Comparison between the OR5 and OR10 cases shows a 17% gain in average rates due to early employment of a variable pressure regime at OR5 (Table 2). At those cases no ice formation was found. At OR5i, only a 4% average rate growth is predicted over a value obtained at OR10i presumably because of extensive ice formation around a wellbore found at the former case. 3.2. Simulations Using Increased Formation Depth. Evolution of ice phase around a wellbore reported at the OR5i case is reminiscent of a secondary hydrate barrier evolution predicted in early simulations.7,8,18 The barrier is initially formed at the top of reservoir and extends down along a wellbore (see discussion about reasons of its appearance in Section 2.5). It has an adverse effect on production rate and causes an increase of gas saturation in a reservoir. Interestingly, in the simulations discussed above in this work no secondary hydrate formation was predicted even for the OR10i and OR5i cases where extensive ice formation was evolved in unit D. In this regard, it is natural to wonder about parameters and conditions of this reservoir model allowing avoid secondary hydrate formation. Can it be attributed to a use of heterogeneous porosity/ saturation and permeability or because of initial conditions of the Mount Elbert reservoir having a temperature profile only a few degrees above the freezing point and, thus, providing limited in situ heat support for the dissociation reaction? Preliminary simulations have shown that an increase of reservoir temperature by one degree can significantly change 1086

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

Figure 10. Production gas rates computed for the OR0, OR0L, OR0i, and OR0Li cases.

production rates41,67 as heat supply to dissociating interface is a major factor affecting production performance.66,68-70 As a consequence, intensive hydrate decomposition promotes faster cooling that might cause hydrate reformation in combination with other factors (Section 2.5). Recently, Moridis et al.67 compared the performance of unit C at the PBU-L106 site in North Slope, Alaska to that of unit D accumulation at the Mount Elbert site during a 2-year production period. They found a superior performance of the former reservoir primarily because of its higher temperature profile exceeding that at the Mount Elbert site by several degrees. To explore an effect of initial conditions on production and possible emergence of gas hydrate reformation, new cases were prepared. Instead of arbitrarily chosen numbers, initial P/T conditions of the reservoir model were set close to those existing at the L-Pad reservoir at Prudhoe Bay Unit on the North Slope Alaska south of the Mount Elbert site. This implies that the reservoir model occurred approximately 180 feet deeper in a warmer temperature range of 3.0-6.5 °C. Assuming the hydrostatic pressure gradient and the local geothermal gradient equal to 0.029 °C/m, the average pressure and temperature were estimated to be 7.40/7.74 MPa and 5.18/6.20 °C for units D and C, respectively. For the layers of the reservoir model reinitialization of initial P/T conditions was performed in response to new formation depth and temperature. In this case, the reservoir exists closer to the base of gas hydrate stability and, thus, can be more easily destabilized by the depressurization method. Two new cases are designated as OR0L and OR0Li, which differ from the OR0 and OR0i ones only by the initial P/T conditions. These new cases also allow us to compare result obtained with permeability based on heterogeneous values and utilizing “default” uniform numbers. Figure 10 displays production gas rates obtained at the OR0, OR0i, OR0L, and OR0Li cases for the first 10 years of production. It can be seen that the reservoir model utilizing increased T and P demonstrates enhanced production performance. Moreover, at ORLi a use of enhanced permeability results in the highest gas rate during first years of production (Figure 10). This prediction is consistent with results of ANOVA and sensitivity analysis of productivity performance showing that interactions between factors responsible for permeability and those controlling heat transfer are the most significant ones at a reservoir scale.68-70 A striking feature of the OR0Li and OR0L production curves depicted in Figure 10 is a cusp appearing at around 4.7 years of production and showing an abrupt decrease in rate. It is

ARTICLE

interesting that at OR0Li, the production rate curve experiences a larger drop than at OR0L. A close inspection of the hydrate saturation distribution revealed secondary hydrate formation evolving around a wellbore at the time that coincides with the production rate drop. Figure 11 presents hydrate and gas saturation distribution around a wellbore for unit D after 5 years of production for both cases. In the OR0Li case, the reformation barrier evolves along a wellbore and effectively shields gas flow. At OR0L, only a smaller part of the wellbore is blocked for gas production and that apparently leads to a higher average production rate. Gas is accumulated behind the barrier and its flow has to divert from the barrier in order to reach a wellbore (Figure 11). Because the only difference between the two cases is a choice of permeability models, that choice controls the appearance of secondary hydrate and, consequently, affects gas production. In other words, a use of heterogeneous intrinsic permeability results in an increase in production rate at “warmer” initial conditions in comparison with the case employing uniform “default” permeability values. It is also clear that a use of warmer and higher pressure initial conditions provides higher production rates that those computed for the shallower (and colder) conditions. However, high production performance is accompanied with evolving secondary hydrate formation around a wellbore, which degrades gas production. These results show the significance of permeability modeling for prediction of hydrate reformation. As this parameter determines an ability of a reservoir to transmit mobile phases, variation of gas saturation and rates along a wellbore could cause different degree of cooling originated from the Joule-Thomson effect. As a result, hydrate reformation phenomenon might not occur at a predetermined region of porous media along a wellbore predicted using the uniform values of permeability (OR0Li) and could follow a different pattern like one observed in the OR0L case. A consequence of that is a more optimistic prediction of production potential. One should bear in mind that predictions of hydrate reformation are made at the equilibrium model. This implies that hydrate (re)formation occurs instantaneously as soon as P/T conditions enter the hydrate stability zone. In reality, kinetics of reformation could be an important factor determining evolution of that process. At the present time, there is a gap in knowledge that hinders the development of a quantitative kinetic model for gas hydrate formation. In part, that is because the hydrate crystallization is a stochastic process exhibiting induction time. That time can vary significantly depending on conditions, water history, chemistry of a hydrate-bearing sediment, etc. Implementation of a model capable to describe kinetics of hydrate (re)formation would be a valuable addition into the current simulator codes. 3.3. Economic Assessment of Productivity Performance. For this assessment, we used an approach similar to one reported by Nyayapathi.70 The net present value (NPV) criterion and the breakeven price of a wellhead gas were computed for cases of interest to make economic evaluations. The breakeven price of gas means that profit becomes zero over a production time period at a chosen discount rate. To calculate a present value at each year over a 25-year production time (it was shown that initial years of returns dominate the net present value of a project70) the following approximate economic data were used: wellhead gas price of $7/mscf; drilling cost of $6 million; fixed cost of $1 million/year; operating and maintenance cost of $0.15/mscf; discount rate of 15%; and corporate tax rate of 35%. These prices are representative of the type of a reservoir 1087

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels

ARTICLE

Figure 11. Secondary hydrate and gas saturation around a wellbore in unit D after 5 years of production at the OR0Li and OR0L cases.

Table 3. Revenue and Breakeven Prices of Wellhead Gas for Cases of Interest

a

OR0

OR10

OR5

OR0i

OR0L

OR0Li

NPV, $M 25 yrs

-10.47

-9.98

-9.44

-8.42

-6.41a

-7.51a

breakeven $/mscf

67

48

36

25

19a

26a

10 Years.

considered here and its geological conditions. The economic viability of a project is subject to its local geological settings and local taxes and tariffs. Table 3 collects net present values and breakeven prices of gas for the OR0/5/10, OR0i, OR0L, and OR0Li cases. The analysis was done for 25 years of production at the OR0/5/10 and OR0i cases using the initial conditions of the Mount Elbert site and for 10 years at the cases utilizing the initial conditions of PBU L-Pad. The data of Table 3 show that gas production from the modeled reservoir remains nonprofitable even with increased productivity owing to lowering bottomhole pressure below the quadruple point at the OR5 and OR10 cases. At the OR0i case, the use of enhanced permeability provides better performance over OR0 and the breakeven price of gas decreases from $67 to $25/mscf. It is noted that similar permeability enhancement resulted in worse performance for OR0Li in comparison with the OR0L case. That is a direct consequence of gas flow blockage due to extended secondary hydrate formation around a wellbore at the case adopting uniform permeability values for the hydratebearing units.

4. CONCLUSION Numerical simulations were used to evaluate gas production potential of the reservoir based on the Mount Elbert methane hydrate deposit at Milne Point on Alaska’s North Slope. The advanced reservoir model constituting hydrate-bearing units D and C separated by an intermediate silt layer was prepared using vertically varying porosity, hydrate saturation, and irreducible water saturation derived from well log data. The permeabilityporosity model was implemented to provide variable permeability values. The depressurization method was used to induce gas hydrate decomposition at a constant bottomhole pressure. To improve production performance the variable bottomhole pressure regimes were successfully employed by slowly reducing pressure from 2.7 to 2.0 MPa over a 15-year period. At the most promising scenario the pressure was reduced below the quadruple point after 5 years of production and the reservoir provided a sustained gas production rate over a 30-year period by using a conventional well design. The predicted total volume of produced gas was nearly doubled in comparison with amount of gas produced at the traditional constant pressure regime. No extensive ice formation was predicted during a 30-year production period. That finding was attributed to a use of heterogeneous permeability values deduced from an experimental permeability-porosity relationship. If those values are replaced with uniform permeability values assigned to hydratebearing units, it results in higher production rates and also promotes formation of a multimeter ice block in unit D. This ice formation is evolved in the middle section of the unit and does 1088

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels not represent a significant obstacle for gas flowing to a wellbore. This shows that a use of permeability values coupled with porosity is important to predict possible ice formation in a producing reservoir. Secondary hydrate formation around a wellbore is known as a phenomenon that can severely hinder gas production. Hydrate reformation was not confirmed during gas production from the reservoir model utilizing cold initial conditions at the Mount Elbert site. Special cases adopting warmer and higher pressure initial conditions revealed enhanced production rates and development of hydrate reformation in the upper unit D. If heterogeneous permeability values are used hydrate reformation occurs only in a localized region that alleviates its adverse effect on gas production. These results emphasize the significance of using heterogeneous porosity, hydrate saturation, and permeability to reliably predict production potential of a reservoir and phenomena affecting production rate. Strong sensitivity of the productivity performance to a choice of permeability values and initial temperature is expected since mass- and heat-transfer are key factors controlling production at a reservoir scale. Economic estimates of the reservoir performance revealed that the reduction of bottomhole pressure below the quadruple point alone cannot bring the production of the relatively cold gas hydrate reservoir into a profitable zone. However, as a tool to enhance production rates at reservoirs showing limited performance and producing at a breakeven price that approach can be a viable solution. It is also possible that pressure reduction could deliver large benefits in terms of production rates at high performance reservoirs and in conjunction with more complex production scenarios and well designs.

’ ASSOCIATED CONTENT

bS

Supporting Information. Porosities, hydrate saturations, irreducible water saturations, and intrinsic permeabilities assigned to layers of hydrate-bearing units D and C; phase diagram of the water-CH4-hydrate system; graphs of relative permeabilities vs water saturation and capillary pressure vs water saturation; hydrate saturation, gas saturation, and pressure distribution for the EVS model after 20 year of production at three different scenarios. This information is available free of charge via the Internet at http://pubs.acs.org

’ AUTHOR INFORMATION Corresponding Author

*Phone: þ1 412 386 5741. Fax: þ1 412 386 5920. E-mail: [email protected].

’ DISCLOSURE Reference in this report to any specific product, process, or service is to facilitate understanding and does not imply its endorsement or favoring by the United States Department of Energy. ’ ACKNOWLEDGMENT E.M. and B.A. performed this work under contract DEFE0004000, Subtask 4000.4.605.261.001 in support of the National Energy Technology Laboratory’s Office of Research and Development. We thank R. Hunter (ASRC Energy) for sharing permeability data used in the paper and W. Winters (U.S.

ARTICLE

Geological Survey) for providing valuable comments for the Simulation details section.

’ NOMENCLATURE kdry = dry thermal conductivity, W/m K kwet = wet thermal conductivity, W/m K kΘ = composite thermal conductivity, W/m K k = intrinsic permeability, m2 kr = relative permeability, m2 P = pressure, Pa Pinit = initial pressure, Pa Pw = pressure at the wellbore, Pa PQ = equilibrium pressure at the quadruple point, Pa ΔPw = pressure drop at the wellbore, Pa SA = aqueous saturation SmaxA = maximum aqueous saturation SG = gas saturation SH = hydrate saturation SI = ice saturation SirrA = irreducible aqueous saturation SirrG = irreducible gas saturation T = temperature, K Greek letters

R = pore compressibility, Pa-1 φ = porosity λ = parameter for the capillary pressure model Subscript and Superscript

A = aqueous c = critical e = equilibrium G = gas H = hydrate I = ice init = initial irr = irreducible max = maximum n = exponent o = free of hydrate and ice conditions p = pressure r = relative w = conditions at the wellbore Θ = composite # = adjusted

’ REFERENCES (1) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gas, 3rd ed.; CRC Press: Boca Raton, FL, 2008. (2) Boswell, R. J. Pet. Sci. Eng. 2007, 56, 9–13. (3) Goel, N. J. Pet. Sci. Eng. 2006, 51, 169–184. (4) White, M. D.; McGrail, B. P. Numerical simulation of methane hydrate production from geologic formations via carbon dioxide injection. Presented at the 2008 Offshore Technology Conference, Houston, TX, May 5-8, 2008; Paper OTC 19458. (5) White, M .D.; Wurstner, S. K.; McGrail, B. P. Mar. Pet. Geol. 2009No. http://dx.doi.org/10.1016/j.marpetgeo.2009.06.008. (6) Moridis, G. J.; Collett, T.; Dallimore, S.; Satoh, T.; Hancock, S.; Weatherhill, B. J. Pet. Sci. Eng. 2004, 43, 219–239. (7) Moridis, G. J.; Kowalsky, M. B.; Pruess., K. Depressurizationinduced gas production from Class 1 hydrate deposits. Presented at the 2005 Society of Petroleum Engineers (SPE) Annual Technical 1089

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels Conference and Exhibition, Dallas, TX, October 9-12, 2005; Paper SPE 97266. (8) Moridis, G. J.; Reagan, M. T. Strategies for gas production from oceanic Class 3 hydrate accumulations. Presented at the 2007 Offshore Technology Conference, Houston, TX, April 30-May 3, 2007; Paper OTC 18865. (9) Moridis, G. J.; Collett, T.; Boswell, R.; Kurihara, M.; Reagan, M. T.; Koh, C.; Sloan, E. D. Toward production from gas hydrates: current status, assessment of resources, and simulation-based evaluation of technology and potential. Presented at the 2008 Society of Petroleum Engineers (SPE) Unconventional Reservoirs Conference, Keystone, CO, February 10-12, 2008; Paper SPE 114163. (10) Kurihara, M.; Sato, A.; Ouchi, H.; Narita, H.; Masuda, Y. Prediction of gas productivity from Eastern Nankai Trough methanehydrate reservoirs. Presented at the 2008 Offshore Technology Conference, Houston, TX, May 5-8, 2008; Paper OTC 19382. (11) Kurihara, M.; Funatsu, K.; Ouchi, H.; Masuda, Y.; Yasuda, M.; Yamamoto, K.; Numasawa, M.; Fujii, T.; Narita, H.; Dallimore, S. R.; Wright, F. Analysis of the JOGMEC/NRCAN/AURORA MALLIK has hydrate production test through numerical simulation. Proceedings of 6th International Conference on Gas Hydrates, Vancouver, British Columbia, Canada, July 6-10, 2008. (12) Kurihara, M.; Sato, A.; Ouchi, H.; Narita, H.; Ebinuma, T.; Suzuki, K.; Masuda, Y. Prediction of production test performances in Eastern Nankai Trough methane hydrate reservoirs using 3D reservoir model. Presented at the 2010 Offshore Technology Conference, Houston, TX, May 3-6, 2010; Paper OTC 20737. (13) Anderson, B. J.; Wilder, J. W.; Kurihara, M.; White, M. D.; Moridis, G. J.; Wilson, S. J.; Pooladi-Darvish, M.; Masuda, Y.; Collett, T. S.; Hunter, R. B.; Narita, H.; Rose, K.; Boswell, R. Analysis of modular dynamic formation test results from the Mount Elbert-01 stratigraphic test well, Milne Point Unit, North Slope, Alaska. Proceedings of 6th International Conference on Gas Hydrates, Vancouver, British Columbia, Canada. July 6-10, 2008. (14) Anderson, B. J. Effects of reservoir heterogeneity on productivity of gas hydrate reservoirs. Fire in the Ice, Methane Hydrate Newsletter, National Energy Technology Laboratory, Fall 2008; http://www. netl.doe.gov/technologies/oil-gas/publications/Hydrates/Newsletter/ HMNewsFall08.pdf. (15) Anderson, B. J.; Kurihara, M.; White, M. D.; Moridis, G. J.; Wilson, S. J.; Pooladi-Darvish, M.; Gaddipati, M.; Masuda, Y.; Collett, T. S.; Hunter, R. B.; Narita, H.; Rose, K.; Boswell, R. Mar. Pet. Geol. 2010No. http://dx.doi.org/10.1016/j.marpetgeo.2010.01.015. (16) Pooladi-Darvish, M.; Hong, H. Mar. Pet. Geol. 2010No. http:// dx.doi.org/10.1016/j.marpetgeo.2010.01.006. (17) Reagan, M. T.; Moridis, G. J.; Zhang, K. Sensitivity analysis of gas production from Class 2 and Class 3 hydrate deposits, Presented at the 2008 Offshore Technology Conference, Houston, TX, May 5-8, 2008; Paper OTC 19554. (18) Moridis, G. J.; Reagan, M. T. Gas production from oceanic Class 2 hydrate accumulations. Presented at the 2007 Offshore Technology Conference, Houston, TX, April 30-May 3, 2007; Paper OTC 18866. (19) Kurihara, M.; Ouchi, H.; Inoue, T.; Yonezawa, T.; Masuda, Y.; Dallimore, S. R.; Collett, T. S. Analysis of the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate thermal-production test through numerical simulation. In Scientific Results from the Mallik 2002 Gas Hydrate Production Research Well Program, Mackenzie Delta, Northwest Territories, Canada; Dallimore, S.R., Collett, T.S., Eds.; Bulletin 585; Geological Survey of Canada, 2005. (20) Alp, D.; Parlaktuna, M.; Moridis, G. J. Energy Convers. Manage. 2007, 48, 1864–1879. (21) White, M. D.; Oostrom, M. STOMP Subsurface Transport Over Multiple Phase: User’s Guide Report; PNNL-15782 (UC-2010); Pacific Northwest National Laboratory: Richland, WA, 2006. (22) MH-21 Research Consortium, 2002; http://www.mh21japan. gr.jp/english. (23) CMG STARS. Computer Modeling Group Ltd, Calgary, Alberta, Canada, 2007; http://www.cmgroup.com/software/stars.htm.

ARTICLE

(24) Wilson, S. J.; Hunter, R. B.; Collett, T. S.; Hancock, S.; Boswell, R.; Anderson, B. J. Mar. Pet. Geol. 2010No. http://dx.doi.org/10.1016/j. marpetgeo.2010.03.007. (25) Moridis, G. J.; Kowalsky, M. B.; Pruess, K. TOUGHþHYDRATE v1.0 User’s Manual: A code for the simulation of system behavior in hydrate-bearing geologic media; Report LBNL-149E; Lawrence Berkeley National Laboratory: Berkeley, CA, 2008; http://esd.lbl.gov/toughþ/ manuals/TplusH_Manual_v1.pdf. (26) National Energy Technology Laboratory. The National Methane Hydrates R&D Program Hydrate Modeling, 2010; http://www.netl.doe.gov/ technologies/oil-gas/FutureSupply/MethaneHydrates/rd-program/Tough FX/ToughFx.html ; http://www.netl.doe.gov/technologies/oil-gas/FutureSupply/MethaneHydrates/MH_CodeCompare/MH_CodeCompare.html. (27) Anderson, B.; Hancock, S.; Wilson, S.; Enger, C.; Collett, T.; Boswell, R.; Hunter, R Mar. Pet. Geol. 2010No. http://dx.doi.org/ 10.1016/j.marpetgeo.2010.02.012. (28) Wilder, J.; Moridis, G. J.; Willson, S.; Kurihara, M.; White, M. D.; Masuda, Y.; Anderson, B. J.; Collette, T. S.; Hunter, R. B.; Narita, H.; Pooladi-Darvish, M.; Rose, K.; Boswell, R. An international effort to compare gas hydrate reservoir simulators. Proceedings of 6th International Conference on Gas Hydrates, Vancouver, British Columbia, July 7-11, 2008. (29) Dallimore, S. R., Collett, T. S., Eds. Scientific Results from the Mallik 2002 Gas Hydrate Production Well Program, Mackenzie Delta, Northwest Territories; Bulletin 585; Geological Survey of Canada, 2005. (30) Dallimore, S.; Wright, J. F.; Nixon, F. M.; Kurhiara, M.; Yamamoto, K.; Fugjii, T.; Fujii, K.; Numasawa, M.; Yasuda, M.; Imasato, Y. Geologic and porous media factors affecting the 2007 production response characteristics of the JOGMEC/NRCAN/Aurora Mallik gas hydrate production research well. Proceedings of 6th International Conference on Gas Hydrates, Vancouver, British Columbia, Canada, July 6-10, 2008. (31) Yamamoto, K.; Dallimore, S. Aurora-JOGMEC-NRCan Mallik 2006-2008 Gas Hydrate Research Project Progress. Fire in the Ice, Methane Hydrate Newsletter, National Energy Technology Laboratory, Summer 2008; http://www.netl.doe.gov/technologies/oil-gas/publications/Hydrates/Newsletter/HMNewsSummer08.pdf. (32) Hunter, R. B.; Collett, T. S.; Boswell, R.; Anderson, B. J.; Digert, S. A.; Pospisil, G.; Baker, R.; Weeks, M. Mar. Pet. Geol. 2010No. http:// dx.doi.org/10.1016/j.marpetgeo.2010.02.015. (33) Reagan, M. T.; Kowalsky, M. B.; Moridis, G. J. The effect of reservoir heterogeneity on gas production from hydrate accumulations in the permafrost. Presented at the 2010 Society of Petroleum Engineers (SPE) Western Regional Meeting, Anaheim, CA, May 27-29, 2010; Paper SPE 97266. (34) Collett, T. S.; Lewis, R. E.; Winters, W. J.; Lee., M. W.; Rose, K. K.; Boswell, R. M. Mar. Pet. Geol. 2010No. http://dx.doi.org/ 10.1016/j.marpetgeo.2010.03.016. (35) Lee, M. W.; Collett, T. S. Mar. Pet. Geol. 2009No. http://dx.doi. org/10.1016/j.marpetgeo.2009.06.007. (36) Boswell, R.; Rose, K.; Collett, T. S.; Lee, M.; Winters, W.; Lewis, K. A.; Agena, W. Mar. Pet. Geol. 2009No. http://dx.doi.org/ 10.1016/j.marpetgeo.2009.12.004. (37) Rose, K.; Boswell, R.; Collett, T. Mar. Pet. Geol. 2010No. http://dx.doi.org/10.1016/j.marpetgeo.2010.02.001. (38) Torres, M. E.; Collett, T. S.; Rose, K. K.; Sample, J. C.; Agena, W. F.; Rosenbaum, E. J. Mar. Pet. Geol. 2009No. http://dx.doi.org/ 10.1016/j.marpetgeo.2009.10.001. (39) Lorenson, T. D.; Collett, T. S.; Hunter, R. B. Mar. Pet. Geol. 2010No. http://dx.doi.org/10.1016/j.marpetgeo.2010.02.007. (40) Winters, W.; Walker, M.; Hunter, R.; Collett, T.; Boswell, R.; Rose, K.; Waite, W.; Torres, M.; Patil, S.; Dandekar, A. Mar. Pet. Geol. 2010No. http://dx.doi.org/10.1016/j.marpetgeo.2010.01.008. (41) Moridis, G. J.; Silpngarmlert, S.; Reagan, M. T.; Collett, T.; Zhang, K. Mar. Pet. Geol. 2010No. http://dx.doi.org/10.1016/j. marpetgeo.2010.01.005. (42) Makogon, Y. F. Hydrates of Hydrocarbons, 1st ed.; Penn Well Publishing Company: Tulsa, OK, 1997. 1090

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091

Energy & Fuels (43) van Genuchten, M. T. Soil Sci. Soc. 1980, 44, 892–900. (44) Stone, H. L. Trans. SPE AIME 1970, 249, 214–220. (45) Moridis, G.; Pruess, K. Flow and Transport Simulations Using T2CG1, a Package of Conjugate Gradient Solvers for the TOUGH2 Family of Codes; Report LBL-36235; Lawrence Berkeley Laboratory: Berkeley, CA, 1995. (46) Moridis, G.; Pruess, K. Geothermics 1998, 27, 415–444. (47) Wright, J. F.; Dalimore, S. R.; Nixon, F. M. Influence of grain size and salinity on pressure-temperature thresholds for methane hydrate stability in JAPEX/JNOC/GSC Mallik 2L-38 gas hydrate research-well sediments. In Scientific Results from JAPEX/JNOC/GSC Mallik 2L-38 Gas Hydrate Research-Well Mackenzie, Delta, Northwest Territories, Canada; Dallimore, S. R., Uchida, T., Collett, T. S., Eds.; Bulletin 544; Geological Survey of Canada, 1999. (48) Collett, T.; Boswell, R. The identification of sites for extendedterm gas hydrate reservoir testing on the Alaska North Slope. Fire in the Ice. Methane Hydrate Newsletter, Summer 2009; http://www.netl.doe. gov/technologies/oil-gas/publications/Hydrates/Newsletter/ MHNewsSummer09.pdf. (49) Dickens, G. R.; Quinby-Hunt, M. S. Methane hydrate stability in pore water: a simple theoretical approach for geophysical applications. J. Geophys. Res. 1997, 102, 773–783. (50) Berg, R. R. Method for determining permeability from reservoir rock properties. Trans. - Gulf Coast Assoc. Geol. Soc. 1970, 20, 303–335. (51) Moridis, G. J. SPE J. 2003, 32, 359–370. (52) Johnson, A.; Patil, S.; Dandekar, A. Mar. Pet. Geol. 2009No. http://dx.doi.org/10.1016/j.marpetgeo.2009.10.013. (53) Verma, A.; Pruess, K. J. Geophys. Res. 1988, 93, 1159–1173. (54) Xu, T.; Ontoy, Y.; Molling, P.; Spycher, N.; Parini, M.; Pruess, K. Geothermics 2004, 33, 477–491. (55) Phillips, O. M. Flow and Reactions in Permeable Rocks; Cambridge University Press: Cambridge, U.K. and New York, 1991. (56) Pape, H.; Clauser, C.; Iffland, J. Geophysics 1999, 64, 1447– 1460. (57) Leverett, M. C. Trans. Soc. Pet. Eng. AIME 1941, 142, 152. (58) Clennell, M. B.; Hovland, M.; Booth, J. S.; Henry, P.; Winters, W. J. J. Geophys. Res. 1999, 104, 22985–23003. (59) Tohidi, B.; Anderson, R.; Clennell, M. B.; Burgass, R. W.; Bidercab, A. B. Geology 2001, 29, 867–870. (60) Kleinberg, R. L.; Flaum, C.; Griffin, D. D.; Brewer, P. G.; Malby, G. E.; Peltzer, E. T.; Yesinowski, G. P. J. Geophys. Res. 2003, 108, 2508–2525. (61) Lee, M. W.; Collett, T. S. Geophysics 2001, 66, 763–777. (62) Lee, M. W. Models for Gas Hydrate-Bearing Sediments Inferred from Hydraulic Permeability and Elastic Velocities; USGS Scientific Investigations Report 2008-5219; U.S. Geological Survey, 2008. (63) Tsimpanogiannis, I. N.; Lichtner, P. C. Phys. Rev. E 2006, 74, 056303. (64) Moridis, G. J.; Reagan, M. T.; Zhang, K. The use of horizontal wells in gas production from hydrate accumulations. Proceedings of 6th International Conference on Gas Hydrates, Vancouver, British Columbia, Canada, July 6-10, 2008. (65) Kowalsky, M.; Moridis, G. Energy Convers. Manage. 2007, 48, 1850–1863. (66) Konno, Y.; Masuda, Y.; Oyama, H.; Kurihara, M.; Ouchi, H. Numerical analysis on the rate-dtermining factors of depressurizationinduced gas production from methane hydrate cores. Presented at the 2010 Offshore Technology Conference, Houston, TX, May 3-6, 2010; Paper OTC 20591. (67) Moridis, G. J.; Reagan, M T.; Boyle, K.; Zhang, K. Evaluation of a deposit in the vicinity of the PBU L-106 Site, North Slope, Alaska, for a potential long-term test of gas production from hydrates. Presented at the 2010 Society of Petroleum Engineers (SPE) Western Regional Meeting, Anaheim, CA, May 27-29, 2010; Paper SPE 133601. (68) Myshakin, E.; Gamwo, I.; Warzinski, R. Experimental design applied to simulation of gas productivity performance at reservoir and laboratory scales utilizing factorial ANOVA methodology. Proceedings of the 2009 Toughþ Symposium, Berkeley, CA, September, 14-16, 2009.

ARTICLE

(69) Gaddipati, M.; Nyayapathi, N.; Anderson, B. J. Sensitivity analysis of methane hydrate reservoirs: effects of reservoir parameters on gas productivity and economics; West Virginia STAR Symposium, April, 2009. (70) Nyayapathi, L. Performance and economics of methane hydrate reservoirs. M.S. Thesis, College of Engineering and Mineral Resources, West Virginia University, 2010.

1091

dx.doi.org/10.1021/ef101407b |Energy Fuels 2011, 25, 1077–1091