Mathematical Modeling and Numerical Simulation of Methane

---

Ind. Eng. Chem. Res. 2010, 49, 5231–5245

5231

Mathematical Modeling and Numerical Simulation of Methane Production in a Hydrate Reservoir Isaac K. Gamwo*,† and Yong Liu†,‡ United States Department of Energy, National Energy Technology Laboratory, Pittsburgh, PennsylVania 15236-0940, and Department of Chemical and Petroleum Engineering, UniVersity of Pittsburgh, Pittsburgh, PennsylVania 15261

Methane hydrate, a potential future energy resource, is known to occur naturally in vast quantities beneath the ocean floor and in permafrost regions. It is important to evaluate how much methane is recoverable from these hydrate reserves. This article introduces the theoretical background of HydrateResSim, the National Energy Technology Laboratory (NETL) methane production simulator for hydrate-containing reservoirs, originally developed for NETL by Lawrence Berkeley National Laboratory (LBNL). It describes the mathematical model that governs the dissociation of methane hydrate by depressurization or thermal stimulation of the system, including the transport of multiple temperature-dependent components in multiple phases through a porous medium. The model equations are obtained by incorporating the multiphase Darcy’s law for gas and liquid into both the mass component balances and the energy conservation equations. Two submodels in HydrateResSim for hydrate dissociation are also considered: a kinetic model and a pure thermodynamic model. Contrary to more traditional reservoir simulations, the set of model unknowns or primary variables in HydrateResSim changes throughout the simulation as a result of the formation or dissociation of ice and hydrate phases during the simulation. The primary variable switch method (PVSM) is used to effectively track these phase changes. The equations are solved by utilizing the implicit time finite-difference method on the grid system, which can properly describe phase appearance or disappearance as well as the boundary conditions. The Newton-Raphson method is used to solve the linear equations after discretization and setup of the Jacobian matrix. We report here the application of HydrateResSim to a three-component, four-phase flow system in order to predict the methane produced from a laboratory-scale reservoir. The first results of HydrateResSim code in a peer-reviewed publication are presented in this article. The numerical solution was verified against the state-of-the art simulator TOUGH+Hydrate. The model was then used to compare two dissociation theories: kinetic and pure equilibrium. Generally, the kinetic model revealed a lower dissociation rate than the equilibrium model. The hydrate dissociation patterns differed significantly when the thermal boundary condition was shifted from adiabatic to constant-temperature. The surface area factor was found to have an important effect on the rate of hydrate dissociation for the kinetic model. The deviation between the kinetic and equilibrium models was found to increase with decreasing surface area factor. 1. Introduction 1

The interest in applying computational fluid dynamics (CFD) for the simulation of nonisothermal multiphase multicomponent processes in the subsurface has significantly increased during the past several decades. For example, in the petroleum industry, CFD models have been developed since the 1970s to help optimize oil production by steam flooding. A good overview of the physical and mathematical aspects for the description of numerical reservoir simulation is given in the book by Peaceman.2 Since the 1980s, an increasing number of problems in environmental engineering, such as the contamination of groundwater due to subsurface leakage of petroleum products, has been a concern for governments and industries that has led to the development of multiphase multicomponent models to simulate the transport of contaminants in the subsurface. More recently, the dissociation and production of methane from methane-hydrate-containing reservoirs has posed an important problem requiring the application of nonisothermal, multiphase, multicomponent models. Large amounts of methane, the main component of natural gas, are trapped as methane hydrate.3,4 * To whom correspondence should be addressed. E-mail: gamwo@ netl.doe.gov. Tel.: 412 386-6537. † National Energy Technology Laboratory. ‡ University of Pittsburgh.

During the mid-1990s, the U.S. Geological Survey released reports estimating the natural gas potential of methane hydrate as close to 320000 trillion cubic feet5 compared to 5 trillion cubic feet that make up the world’s currently known gas reserve. The combustion of methane also yields 44% less CO2 than coal, per unit energy release, and 29% less than oil.6 In hydrate, the guest gas molecules are trapped in the cavity formed by the water molecules linked through hydrogen bonding.4 Gas hydrates are stable under relatively high pressures and low temperatures and are found in the natural state below permafrost such as the Mallik gas hydrate field in Canada4 or under the sea floor near the continental margins such as in the Gulf of Mexico.4 Methane hydrate dissociation and formation can be represented by the following reaction4 CH4 + NHH2O(solid) a CH4(gas) + NHH2O(water or ice) (1) where NH, the hydrate number, is the number of water molecules per methane molecule (usually 6.0) in the hydrate lattice. Several conventional methods and their combinations have been proposed for the production of natural gas from hydratecontaining reservoirs.4 The first and most popular method is depressurization, in which the pressure in the reservoir is

10.1021/ie901452v  2010 American Chemical Society Published on Web 03/10/2010

5232

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

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 is inhibitor injection, in which methanol or brine is typically used to shift the thermodynamic equilibrium and dissociate the hydrate. It is also thermodynamically possible to produce methane from hydrate by displacing it with carbon dioxide.7 The exothermic heat released by carbon dioxide hydrate formation is used to overcome the endothermic heat needed for the methane hydrate dissociation. This fourth, unconventional method also offers the potential for carbon sequestration. Over the past two decades or so, several mathematical models and numerical codes have been applied for methane hydrate decomposition and methane production.8-35 Reservoir simulator researchers are in need of an open-source, efficient multiphase multicomponent software to simulate such complex physical and chemical processes. Among the numerical codes, HydrateResSim17 (HRS) is the sole open-source code available to the public (www.netl.doe.gov) that describes hydrate decomposition and release in a hydrate reservoir. HydrateResSim was developed for the National Energy Technology Laboratory (NETL) in 2006 by the Lawrence Berkeley National Laboratory (LBNL). However, numerous undetected bugs in the code and a lack of detailed code-related documents made it extremely difficult to use. To illustrate, over 70 research institutions worldwide have requested and received HRS from NETL; however, to our knowledge, there are no peerreviewed publications to date based on its use, likely owing to the difficulties in using the code and the lack of documentation. Here, we describe the first detailed use of this improved code to model methane production from laboratory-scale hydrate reservoirs under two sets of scenarios, namely, adiabatic and constant-temperature boundary conditions. Both kinetic and equilibrium models were used to simulate methane hydrate dissociation. The mathematical model used in the HydrateResSim code is described in detail, including the development of the governing equations, constitutive equations, boundary and initial conditions, and numerical techniques. A parametric study shows that the rate of hydrate dissociation depends strongly on the surface area factor. 2. Mathematical Modeling For a detailed description of hydrate dissociation, the flow of liquid and gas, and the heat transfer in a multiphase, multicomponent subsurface system, it is necessary to formulate for each phase the mass balance equation, the momentum equation, and the energy balances, as well as the mass balance for each component. However, such a full multiphase, sophisticated model concept, which can effectively reflect complexities in large subsurface systems, will result in excessively large systems of highly coupled partial differential equations with large numbers of unknowns and uncertainties that are difficult to quantify. In the use of HydrateResSim, several simplifications are often made just for practical purposes. For example, a macroscopic description of slow flow processes in porous media, Darcy’s law in its modified form for multiphase flow, is commonly used for gas and liquid flow instead of the classical phase momentum equations. This allows for a decoupled calculation of the fluid velocity that can be inserted into the mass balance equation in order to obtain the general multiphase differential equation. Similarly, for a macroscopic description in a representative elementary volume of a porous medium under consideration,

we state that the assumption of local thermal equilibrium is valid because flow velocities are generally small in the unfractured porous medium. Such an assumption allows for the formulation of a single balance of energy equation for the whole fluid-filled porous medium. Thus, the heat-transfer terms between the fluid phases do not apply. We treat the amounts of thermal energy stored in the fluid phases and in the soil grain additively. We further take into account heat transport due to binary diffusion in the gas phase and heat sources or sinks. For the initial model, the phase balances have been ignored for the sake of simplification. The solids are assumed to be immobile. In view of these assumptions, the nonisothermal four-phase (aqueous, A; gaseous, G; hydrate, H; ice, I) three-component (methane, m; water, w; hydrate, h) hydrate system is described here with the mass balances for each component, the Darcy’s laws for aqueous and gaseous phases, and one energy equation. The number of phases existing in a given control volume is not necessarily constant, as the mass transfer between the phases such as ice formation or melting can cause the appearance and disappearance of phases. The HydrateResSim numerical simulator can describe hydrate dissociation with either an equilibrium or a kinetic submodel. We have chosen two methods for hydrate dissociation in this work to illustrate the difference between the submodels. In the following section, the details of these submodels and the numerical techniques are described. The details of the discretization of governing equations for kinetic model, the numerical techniques of the Newton-Raphson method, and the primary switch method for the equilibrium model are presented in the Appendix. 2.1. Kinetic and Equilibrium Submodels. In the kinetic model, the rate of hydrate dissociation or formation is proportional to the driving force, defined as the difference between the gas fugacity and the equilibrium hydrate fugacity. The Kim-Bishnoi equation23,36 expresses the molar rate of gas generated by hydrate dissociation, n˙, as

(

n˙mole ) k0d exp -

∆E f(AHS, φ, SH)[feH(T) - fG] RT

)

(2)

where AHS is the total surface of hydrate, feH is the equilibrium gas-hydrate-liquid fugacity that can be approximated by the equilibrium pressure PeH,11,28 fG is the gas fugacity that is approximated by the gas pressure PG in pore, ∆E is the activation energy, R is the universal gas constant, T is the temperature, and k0d is an intrinsic dissociation constant independent of pressure and temperature. The pure equilibrium model uses the prevalent temperature and pressure in the reservoir along with a look-up table composed of temperature and pressure at the gas-hydrate-liquid equilibrium. If the reservoir pressure or temperature falls below the equilibrium value, then the hydrate dissociates spontaneously. In the equilibrium model, the hydrate is simply a phase, not a component as in the kinetic model. The equilibrium model is inherently generic and relatively easy to implement with faster convergence. The kinetic rate model, on the other hand, is accurate and detailed but is computationally more intensive. It has been reported that both the kinetic and equilibrium submodels in TOUGH+HYDRATE predict comparable results for reservoir-scale simulations; however, for short-term processes at smaller scales, kinetic limitations can be important.32 In using the kinetic submodel, it is assumed that the gas-liquid-solid systems are composed of four phases, aqueous (A), gaseous (G), hydrate (H), and ice (I). These phases are

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

5233

In eq 3, b FA and b FG are the Darcy’s velocities for the aqueous and gaseous phases, respectively. The Darcy’s velocities of the aqueous phase or gaseous phase are averaged over the cross section of the porous medium. For fluid flow in porous media, the gas and aqueous velocities are assumed to obey the multiphase Darcy’s law Figure 1. Schematic of phases and components in the kinetic model.

composed of the components water (w), methane (m), and hydrate (h) as shown in Figure 1. The mass-transfer processes between phases taken into account in the model concept include dissolution, degassing, evaporation, and condensation. Adsorption is neglected. The hydrate and ice phases can appear or disappear during the production of methane hydrate depending on the local temperature and pressure conditions. In the equilibrium model, the rates for ice or hydrate formation or dissociation are assumed to be very high. At a given temperature, if the local pressure drops below the equilibrium pressure, PeH, hydrate is assumed to dissociate into gas and water instantaneously in the equilibrium model. In the kinetic model, the rate for hydrate formation or dissociation is assumed to be proportional to the difference between the equilibrium fugacity, feH, and the local gas fugacity, fG. In the kinetic model, ice melting or formation is still considered as occurring at chemical equilibrium. In the equations herein, all components are written in lower case, whereas all phases are written in upper case, as shown in Figure 1. The focus of this work is to explain the governing equations and the primary switch method; thus, the use of an inhibitor was neglected. The governing equations need to be modified if any inhibitor such as salt is used. Aside from the possible component distributions in Figure 1, the solid phases of sand, ice, and hydrate are assumed to be stationary; therefore, no momentum equation is needed for the solid phases. The gas phase has a water component (w) because of the vaporization of water, and the aqueous phase has a methane component (m) because of the dissolution of methane at high pressure. The physical meanings of symbols are described in the Nomenclature section. 2.2. Governing Equations (Kinetic Model). In the following mass balances, the transport of components can occur by advection with the flowing phases and by diffusion of the components in the gas phase. The term q, which is used to model the sources or sinks of the respective components inside the model domain or at the boundaries, is also taken into account. Hydrate dissociation is considered as a chemical reaction of molar rate n˙. The adsorption of fluid components at the solid matrix and the effects of hysteresis in the constitutive relationships are neglected. In multiphase flow, a phase can be defined as an identifiable class of material that has a particular inertial response to and interaction with the flow and the potential field in which it is immersed. In the kinetic dissociation model, there are four possible phases, namely, the aqueous phase, the gaseous phase, the hydrate phase, and the ice phase, and three components, namely, methane (m), water (w), and hydrate (h). The continuity equation for component methane (m) can be written as

b FA ) - Κabs

(

ΚrAFA (∇PA - FAb g) µA

)

b ΚrGFG b FG ) - Κabs 1 + (∇PG - FGb g) PG µG

(4)

(5)

In eq 5, b is the Klinkenberg constant accounting for gas slippage effects. The continuity equation for the water (w) component can be written as

The continuity equation for the hydrate (h) component is

For the energy balance equation, the assumption of local thermal equilibrium is valid because flow velocities are generally small in porous media. This assumption permits the formulation of a single energy balance equation for the whole fluid-filled porous medium. Dissipation work is neglected because the velocity gradients in porous media flow are rather small. The energy balance equation can be written based on the summation of the enthalpy over all phases as

In eq 8, HR is the rock enthalpy; HA, HG, and HH are given by m HA ) XAwHAw + XAm(HAm + Hsol )

(9)

HG ) XGwHGw + XGmHGm + Hdeparture

(10)

HH ) XHwHHw + XHmHHm + ∆HH0

(11)

and ∆HH0 is the latent heat of hydrate dissociation. HI ) HIw + ∆HI0

(12)

where ∆HI0 is the latent heat of ice dissociation. The heat flux terms include conductive heat transfer, convective heat transfer, and radiative heat transfer as follows b Fθ ) -Kave∇T + fσσ0∇T4 +

∑

Hβb Fβ

β)A,G

All symbols are explained in the Nomenclature section.

where the superscript θ stands for heat in the derivation.

(13)

5234

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

Table 1. Primary Variables (PVs) in the Kinetic Model A PV1 PV2 PV3 PV4

A+G A+H A+H+G A+I+G I+H+A+G

PA XAm SH T

PG SA SH T

PA SA XAm T

PG SA SG T

PG SA SH SG

XGw + XGm ) 1

(14)

In HydrateResSim, several expressions are available for the average heat conductivity. We have four governing equations, eqs 3, 6, 7, and 8 that relate the variable values at different times. Thus, we can specify four primary variables to describe the system at different times, as shown in Table 1. Unlike traditional simulation work, where the unknowns or primary variables remain unchanged, in this case, the primary variables change throughout the simulation depending on the phases present, as described in Table 1. 2.3. Choosing the Primary Variables. In the simulations reported here, all of the variables are known at time t, and the goal is to calculate these variables at the next time t + ∆t. There are four partial differential equations, eqs 3, 6, 7, and 8, that relate variables at different times. Therefore, four primary variables can be selected to be solved. The primary variables are chosen such that other variables can be expressed as functions of the primary variables. For example, assuming that, for a given element i at time t, only the aqueous phase and the hydrate phase exist, the aqueous pressure PA, aqueous-phase saturation SA, mass fraction of methane in the aqueous phase XAm, and temperature T can be chosen as the primary variables to be solved; alternatively, the hydrate saturation SH could be selected instead of the aqueous saturation SA as a primary variable. The effect of choosing different primary variables on the simulation results is discussed in the Appendix in the Newton-Raphson Method section. Table 1 lists the primary variables used in the numerical simulator HydrateResSim for different phase combinations in this work using the kinetic model. In Table 1, for example, when two phases coexist, say, the aqueous and gaseous phases (A + G) as shown in the second column, the gas pressure, PG; the aqueous-phase saturation, SA; the hydrate saturation, SH; and the temperature, T, are chosen as the primary variables. When three phases coexist, say, the aqueous, gaseous, and hydrate phases (A + H + G), the gas pressure, PG; aqueous-phase saturation, SA; gas-phase saturation, SG; and the temperature, T, are chosen as the primary variables. 2.4. Constitutive Relationships. The governing equations presented in the previous section are not sufficient to fully describe a multiphase system. They must be supplemented with equations that describe the constitutive behavior of the individual phases. The components are denoted by k, and the phases are denoted by β in the following equations. The molecular fraction of component k in phase β, xβk , has the following relationship with the mass fraction of component k in phase β, Xβk Xβk

)

k xβk Mweight

∑

k xβk Mweight

(15)

where Mweight is the molecular weight of component k. The sum of the mass fractions of the different components in the aqueous phase is 1; that is

(17)

XHw ) 0.871 and XHm ) 0.129 when the hydrate number, NH, equals 6. XIw ) 1 because the ice phase is composed of water only. The porosity, φ, of the porous medium is defined as φ)

Vvoid Vtotal

(18)

where Vvoid is the void volume contained in the porous medium sample and Vtotal is the total volume of the porous medium sample. In the numerical simulator HydrateResSim, the porosity of the sand is assumed to be a function of the pressure and temperature; that is, the compressibility and expansivity are considered φ ) fφ(T, P)

(19)

When the void space of the porous medium is filled by two or more immiscible aqueous or gaseous fluids, the saturation at a point with respect to a particular phase is defined as the fraction of the void volume of the porous medium occupied by that particular phase within a representative elementary volume (REV) around the considered point as shown in the equation Sβ )

Vβ Vvoid

(20)

where Vβ is the volume of phase β within a REV. The sum of saturations equals 1; that is SA + SG + SH + SI ) 1

(21)

The aqueous-phase relative permeability, ΚrA, is KrA ) frA(SA, SG)

(22)

The gas-phase relative permeability, ΚrG, is KrG ) frG(SA, SG)

(23)

The pressure difference between the gas phase and the aqueous phase is defined as the capillary pressure Pcapillary ) PG - PA ) fcap(SA, SG)

(24)

Because of the complexity of porous media, several relative permeability and capillary functions are available in HydrateResSim. The density of the aqueous phase, FA, is FA ) fdenA(T, PA, XAm, XAw)

k)w,m k

(16)

Similarly, the sum of the mass fractions of the different components in the gas phase is also 1

PG SA SG SI

The average heat conductivity is defined as Kave ) fave(φ, SA, KA, SG, KG, SH, KH, SI, KI, KR)

XAw + XAm ) 1

(25)

and the density of the gas phase, FG, is FG ) fdenG(T, PG, XGm, XGw)

(26)

where fdenA and fdenG are the functions used to calculate the densities of the aqueous and gaseous phases, respectively.

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

For the aqueous phase, we assume that the concentration of the dissolved component methane is low. Thus, dissolution can be described by Henry’s law. Hence, the mass fraction of methane in the aqueous phase is calculated as follows: The Henry’s constant is a function of temperature Hm ) fhenry(T)

(27)

The molecular fraction of methane in the aqueous phase is defined as

xAm )

{

PG if PG e PeH Hm PeH if PG > PeH Hm

(28)

After calculation of the mole fraction of methane in the aqueous phase using eq 28, the mass fraction of methane in the aqueous phase is calculated using eq 15. The equilibrium pressure, PeH, for the hydrate is a function of temperature PeH ) feH(T)

(29)

The molecular fraction of water in the gas phase is defined using the water vapor pressure and the partial pressure rule xGw )

Pvapor PG

Table 2. Reactor Properties, Initial Conditions, and Pertinent Model Parameters parameter height length initial pressure boundary pressure initial temperature initial water saturation in core initial gas saturation in core intrinsic permeability, k specific heat capacity porosity, φ sand thermal conductivity relative permeability n SirA SirG evolving porous media (EPM) model17 capillary pressure model SirA λ SmxA

value 0.0508 m 0.3 m 3.75 MPa 2.84 MPa 2.3 °C 0.351 0.206 0.0967 × 10-12 m2 800 J/(kg °C) 0.182 8.8 W/(m K) n n krA ) (S * A ) , krG ) (S * G) * S A ) (SA - SirA)/(1 - SirA), S* G ) (SG - SirG)/(1 - SirG) 3.0 0.15 0.05 Pcap ) -P0[(S*)-1/λ - 1]-λ S* ) (SA - SirA)/(SmxA - SirA) 0.14 0.46 1.0

saturation of 0.206. Other reactor properties, initial conditions, and pertinent model parameters are listed in Table 2. 3. Results and Discussion

(30)

The vapor pressure of water at a given temperature is calculated as Pvapor ) fvapor(T)

5235

(31)

where fvapor is the function used to calculate the vapor pressure. After the mole fraction of water in the gas phase is calculated using eq 30, then the mass fraction of water in the gas phase can be calculated using eq 15. 2.5. Boundary and Initial Conditions. We tested the two boundary conditions implemented in the numerical code HydrateResSim. The first is the Dirichlet boundary condition.17 This kind of boundary condition keeps the absolute value of a primary variable fixed. For example, the boundary condition for temperature is kept at 2.3 °C, and the pressure at the outlet remains constant at 2.84 MPa. The assignment of initial conditions and initial phase states must correspond to the Dirichlet conditions at the respective boundaries. The remaining boundary has no flow properties. The second is the Neumann boundary condition.17 A Neumann condition allows an assignment of molar or heat fluxes across the respective boundaries. The molar fluxes correspond to the source/sink terms in the molar mass balance equations of the methane and water components. In the simulations reported here, the simulation domain is surrounded by walls; therefore, no mass flux is applied at all the walls. Neumann conditions can also be applied for the energy balance of eq 8. Adiabatic conditions are assumed, and the heat flux is set to zero around the boundaries. This restriction to two kinds of boundary conditions causes difficulties in some practical flow problems involving natural heat convection (cooling law) at the boundaries. Initially, the reservoir is at a pressure of 3.75 MPa, a temperature of 2.3 °C, a hydrate saturation of 0.443, and a gas

Four sets of numerical simulations were performed for the production of methane from a laboratory-scale hydrate reservoir. These simulations were designed with the following four objectives: (1) compare the HydrateResSim kinetic model and equilibrium model results, (2) compare and explain the effects of the adiabatic and constant-temperature boundary conditions on the dissociation patterns, (3) validate the HydrateResSim CFD modeling results with those from the existing state-ofthe-art software TOUGH+HYDRATE for the production of methane from the laboratory-scale hydrate reservoir, and (4) show the sensitivity of the production rate to the hydrate surface area factor. Cumulative methane produced at the simulated well, saturation of all of the phases, and pressure and temperature distributions in the reservoir were predicted for the various production scenarios. All of the simulations were based on a laboratory-scale experiment.24 Figure 2 shows a two-dimensional cross section of the computational domain for the laboratory-scale reservoir, which is 5.08 cm in diameter and 30 cm in length. The domain for the results shown consists of a uniform 100 × 20 grid, with ∆x ) 0.30 cm and ∆y ) 0.25 cm. The reactor is initially filled with methane hydrate at 3.75 MPa and 2.3 °C and with initial saturations of 0.443, 0.351, and 0.206 for the hydrate, aqueous, and gas phases, respectively. The simulated well (outlet) is located at the left center as indicated in red in Figure 2 and has an opening of 0.508 cm. The depressurization technique was used by imposing a pressure of 2.84 MPa and a constant temperature of 2.3 °C at the well, which initiates the dissociation and release of methane. Two thermal boundary conditions were tested at the walls: zero heat flux around the boundaries (adiabatic Neumann conditions) and a constant temperature of 2.3 °C around the boundaries (Dirichlet thermal boundary conditions). No mass flow was allowed across the boundaries except at the well. The initial time step was set at 1 s, and HydrateResSim automatically adjusted the time step based on preset convergence criteria. We monitored the results at three

5236

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

Figure 2. Schematic diagram of the simulated 2D laboratory-scale reservoir (converged grid: 2000 elements).

Figure 3. Comparison of pressure as a function of time between the equilibrium and kinetic models at different observation points for adiabatic boundary conditions.

observations points (shown in Figure 2) during all simulations. The points are located at the height y ) 1.7 cm and the lengths x ) 7.5, 15, and 22.5 cm. The computational settings are summarized in Table 2. Prior to conducting extensive simulations, a search for a mesh-independent solution was conducted using the following meshes: 20 × 20, 50 × 20, 100 × 20, and 50 × 40. The results of the grid-independence study will be reported elsewhere to avoid extending this article. The mesh size was varied in both the x and y directions until the simulations reached a numerical asymptotic solution with 2000 grid elements (100 × 20). Subsequent simulations were conducted with 100 × 20 grid points. The simulation based on the kinetic model23 was carried out and compared with the equilibrium model based on the phase equilibrium curve.32 Figure 3 compares the temporal distribution of pressure at the three observation points using the kinetic and equilibrium models. Both models show that, after the onset of depressurization at the well, the pressure inside the reservoir rapidly drops to the depressurized pressure of 2.84 MPa, inducing the flow of the aqueous and gaseous phases toward the well. However, the equilibrium model pressure drops faster than the kinetic model pressure because we assume that the hydrate dissociates instantaneously in the equilibrium model. Such detailed results are not predicted using the equilibrium model. Figure 4a compares the hydrate saturation distribution using both the HydrateResSim equilibrium and kinetic models with the adiabatic boundary conditions. As the well pressure dropped to 2.84 MPa, the hydrate closest to the well dissociated first; then the dissociation expanded from the well throughout the reservoir. Figure 4a shows that, as soon as dissociation begins, curvature develops in the dissociation front owing to the well opening location at the center of the left face (x ) 0, y ) 0.0254

m as shown in Figure 2). After 66 min, the dissociation front becomes more uniform in the y direction and continues to progress in this manner for the rest of the simulation. Figure 4a shows that the two models predict the same behavior; however, the transition of hydrate saturation at the dissociation front is sharper in the equilibrium model than in the kinetic model. This phenomenon is explained by the fact that hydrate dissociates instantaneously in the equilibrium model whereas it dissociates gradually in the kinetic model. Figure 4b shows a comparison of the hydrate saturation as a function of time between the equilibrium and kinetic models at the three observation points. The drop in hydrate saturation is more abrupt in the equilibrium model compared to that observed with the kinetic model. Assuming that the dissociation front propagation is linear, the propagation velocity of this moving front was calculated to be approximately 3 µm s-1. Figure 4c presents a comparison of simulations performed using HydrateResSim for both the equilibrium and kinetic models with constant-temperature boundary conditions. It also shows a comparison of the kinetic results with those obtained using the well-known hydrate reservoir simulator TOUGH+ HYDRATE19 that was also developed at Lawrence Berkeley National Laboratory.37 The constant temperature (2.3 °C) at the boundaries permits continuous heat transfer to the simulation domain as hydrate dissociates. In Figure 4c, in contrast to the previous adiabatic case in Figure 4a, the hydrate close to the boundaries dissociates faster than the hydrate inside the domain owing to the heat transfer from the boundary. Figure 4c shows that, when the equilibrium model is used, the hydrate dissociates similarly to melting ice, shrinking from all directions, whereas it dissociates throughout the whole reservoir with the kinetic model. As previously mentioned, the hydrate is assumed to dissociate instantaneously in the equilibrium model, whereas the hydrate dissociates gradually in the kinetic model. In actual laboratory-scale experiments, hydrate dissociation is normally assumed to be kinetically limited because of the heat available from the surroundings (e.g., the size of the experimental core holder is small compared with the size of experimental room), whereas heat transfer is considered as the limiting process for gas hydrate dissociation in field-scale reservoirs. Rather than basing the choice of model solely on the scale of the process, it might be prudent, based on the results in Figure 4a,c, to choose a model based on the anticipated amount of heat transfer. A simple moving-front model might not be suitable when there is enough heat available at the boundaries to promote hydrate dissociation. Figure 4c also shows that HydrateResSim produces results that are very similar to those obtained using the stateof-the-art TOUGH+HYDRATE simulation code even when the kinetic models were used. HydrateResSim therefore appears to be a suitable option for performing simulation research on hydrate-containing systems. Figure 5a shows a comparison of gas saturation distributions between the equilibrium and kinetic models obtained using HydrateResSim at different times under adiabatic boundary

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

5237

Figure 4. (a) Comparison of hydrate saturation distributions between the equilibrium and kinetic models at different times for adiabatic boundary conditions using HydrateResSim. (b) Comparison of hydrate saturation as a function of time between the equilibrium and kinetic models at different observation points for the adiabatic boundary conditions using HydrateResSim. (c) Comparison of hydrate saturation distribution at different times using HydrateResSim and TOUGH+HYDRATE for a constant temperature (2.3 °C) at the boundaries.

5238

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

Figure 5. (a) Comparison of gas saturation distributions between the equilibrium and kinetic models at different times under adiabatic boundary conditions. (b) Comparison of gas saturation as a function of time between the equilibrium and kinetic models at different observation points under adiabatic boundary conditions.

conditions. It shows that the gas saturation has a maximum value at the moving front. The gas saturation decreases from the moving front to the well and is due to the flow of the gaseous phases out of the reservoir. Figure 5b shows a comparison of gas saturation as a function of time between the equilibrium and kinetic models at three observation points. Figure 5b shows that, at a given location, such as x ) 0.15, the gas saturation can be roughly divided into four stages for both the equilibrium and kinetic models. The transition of gas saturation in the kinetic model is smoother than in the equilibrium model. Stage 1 runs from the initial time to 580 min. In this stage, the gas saturation drops continuously because of the flow of gas out of the domain. Stage 2 runs from 580 to 720 min and represents an increase in gas saturation because of hydrate dissociation. Stage 3 runs from 720 to 2400 min. In this stage, the gas saturation is composed of two processes: one is the decrease of gas saturation due to the gas flow to the well opening, and the other is the increase of gas because of the hydrate that is dissociating between x ) 0.15 m and x ) 0.3 m and flowing into the zone at x ) 0.15 m. The combination of these two processes results in the slow or gradual decrease of gas saturation at x ) 0.15 in stage 3. Stage

4 runs from 2400 to 5000 min and represents the completion of hydrate dissociation and the drop in gas saturation as the remainder of the released gas moves to the well opening. Figure 6a shows a comparison of temperature distributions between the equilibrium and kinetic models at different times under adiabatic boundary conditions. For these conditions, parts b and c of Figure 6 show comparisons of temperature as functions of long and short times, respectively, at the three observation points. Prior to the start of dissociation, the temperature of the reservoir is uniform at 2.3 °C. The dissociation is initiated by dropping the pressure to 2.84 MPa. As previously noted, this pressure drop is quickly transmitted throughout the reservoir, which results in hydrate dissociation throughout the reservoir. The endothermic nature of hydrate dissociation rapidly cools the system. Figure 6a shows that the temperature has dropped to 1 °C over most of the domain after 33 min. Only the area near the well opening remains at a higher temperature, due, in part, to heat transfer as the released gas and water flow toward the well opening. However, the temperature at the well opening is artificially maintained at 2.3 °C throughout the simulations and therefore also serves as a point

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

5239

Figure 6. (a) Comparison of temperature distributions between the equilibrium and kinetic models under adiabatic boundary conditions. (b) Comparison of temperature as a function of time between the equilibrium and kinetic models under adiabatic boundary conditions at different observation points. (c) Comparison of temperature as a function of time between the equilibrium and kinetic models under adiabatic boundary conditions at different observation points early in the dissociation process.

5240

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

source of heat that gradually warms the content of the reservoir. To further illustrate the impact of the artificially heated well opening, parts b and c of Figure 6 show the complete and initial temperature histories, respectively, inside the reservoir during the dissociation process. From this information, the temperature as a function of time for the different observation points can be divided into four stages. Stage 1 represents the rapid temperature drop discussed above as hydrate initially dissociates throughout the reactor in response to the pressure drop that initiates the dissociation process. This stage takes about 7 min. In the second stage, the heat from the well opening gradually extends through the reactor, as evidenced by the constant temperature at the three observation points. In this stage, the temperature is relatively constant, as further hydrate dissociation is prevented under the constant-pressure and adiabatic scenario. As the heat from the constant-temperature well opening reaches each observation point, the temperature in the system begins to gradually rise. In this third stage, additional hydrate decomposes, and the shape of the temperature rise profile is a combination of the heat input and the endothermic nature of the hydrate reaction induced by the temperature rise. At about 2400 min, the hydrate dissociation is complete throughout the reactor, and the temperatures gradually rise in the fourth and final stage, as the entire hydratedepleted reservoir comes to equilibrium with the well opening at 2.3 °C. A key point from these observations is that the parameters set for the well opening can dramatically influence the behavior of the system. Particularly, in the absence of heat flow from the surroundings, the heat flow from the well opening can result in the hydrate dissociation process after the initial dissociation period. Figure 7a shows the results of a parametric study of the effect of the surface area factor used in the kinetic model on the cumulative methane produced at the outlet. With decreasing the area factor, the hydrate dissociates at a lower rate and needs longer time to reach steady state. Figure 7b shows a comparison of hydrate saturation at one observation point (x ) 15 cm, y ) 1.7 cm) between the equilibrium model and the kinetic model with different area factors. Figure 7b shows that the hydrate saturation drops spontaneously for the equilibrium model, whereas the hydrate saturation drops gradually for the kinetic model. With a smaller area factor, the hydrate takes longer time to dissociate. Parts a and b of Figure 7 show that the area factor is important to the rate of hydrate dissociation. We suggest more experimental work to measure the area factor correctly. 4. Conclusions We have presented the development of HydrateResSim, a multicomponent, multiphase fluid- and heat-flow simulation code to describe methane production from dissociating hydrate in laboratory-scale systems designed to mimic larger-scale reservoir conditions. This simulation code includes both a kinetic model and a pure equilibrium model for simulation of the hydrate dissociation process. The simulation code has been applied to a three-component, four-phase flow system to predict the methane produced from a laboratory-scale reservoir. These are the first detailed results of HydrateResSim presented in the open, peer-reviewed literature. The major computational findings from the present study are as follows: We showed the details of the primary switch method used in the TOUGH family code for hydrate dissociation and formation. The HydrateResSim results obtained using isothermal boundary conditions were found to be in good agreement withresultsfromthestate-of-the-artmodel,TOUGH+HYDRATE,

Figure 7. (a) Comparison of cumulative gas produced between the equilibrium and kinetic models with different area factors for adiabatic boundary conditions. (b) Comparison of hydrate saturation at (x ) 15 cm, y ) 1.7 cm) between the equilibrium and kinetic models with different area factors for adiabatic boundary conditions.

a fact indicating the validity of this unique freeware, opensource code that is capable of describing methane hydrate behavior. Microscopic examination of the pressure and temperature distributions in the reservoir revealed that the equilibrium model generally overpredicts the rate of hydrate dissociation compared to the kinetic model. The effect of thermal boundary conditions on the dissociation pattern was also clearly demonstrated. Adiabatic boundary conditions lead to a single dissociation front, whereas constanttemperature boundary conditions result in dissociation propagation from all boundaries. Finally, the effect of the artificially maintained temperature in the simulated well opening was clearly shown when adiabatic boundary conditions were employed. This is one of the major challenges in utilizing these simulators to model actual laboratory-scale experimental data. Even in the case where heat transfer is permitted, the simulation conditions and those that exist in the laboratory might be significantly different to permit accurate simulations of experimental results. Perhaps this latter case can be resolved by more complicated simulation models that better represent laboratory conditions and laboratory-reactor-scale configurations. We are endeavoring to make such improvements. A parametric study showed that the hydrate dissociation rate decreases with

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

increasing surface area factor. We strongly recommend experimental investigations to measure the surface area factor accurately. Future work is planned to apply the model to experimental and field data. Acknowledgment The authors thank Professor Hamid Arastoopour of Illinois Institute of Technology and Professor Goodarz Ahmadi of Clarkson University for their comments and suggestions. This research was supported in part by an appointment to the U.S. Department of Energy (DOE) Postgraduate Research Program at the National Energy Technology Laboratory administrated by the Oak Ridge Institute for Science and Education. We thank Ray Boswell for recommending HydrateResSim and Jamie Brown for her continuous support throughout the project. We thank Robert P. Warzinski of the National Energy Technology Laboratory (NETL) for his candid comments. We also thank Evgeniy M. Myshakin of NETL and George J. Moridis of Lawrence Berkeley National Laboratory for their assistance during this project.

w ) water component θ ) energy component Subscripts A ) aqueous phase cap ) capillary pressure G ) gaseous phase H ) hydrate phase I ) ice phase nm ) element n with direction m R ) rock phase rβ ) relative permeability for phase β β ) phase β (β could be gas, aqueous, hydrate, ice, or sand)

Discretization of Governing Equations for the Kinetic Model The governing eqs 3, 6, 7, and 8 can be rewritten in a compact form to facilitate the discretization as follows:29 The continuity equation for the methane component, k ) m, is

Nomenclature Anm ) surface area of element n with direction m (m2) f ) function Fk ) flux of component k fσ ) radiance emittance factor g ) gravitational constant Hβ ) enthalpy of phase β (J mol-1 or J kg-1) Jβk ) diffusion of component k in phase β Κabs ) absolute permeability (m2) k0d ) intrinsic dissociation constant (mol m-2 Pa-1 s-1) ΚrA ) relative permeability for the aqueous phase ΚrG ) relative permeability for the gaseous phase Mk ) mass of component k (kg) k Mweight ) molecular weight of molecular k (kg mol-1) PeH ) equilibrium pressure of the hydrate phase (Pa) PeI ) equilibrium pressure of the ice phase (Pa) Pβ ) pressure of phase β (Pa) qk ) source or sink of component k due to the flow of gas or aqueous phase Sβ ) saturation of phase β t ) time (s) xβk ) mole fraction of component k in phase β Xβk ) mass fraction of component k in phase β

Superscripts h ) hydrate component k ) component k (k could be methane, water, hydrate, or energy) m ) methane component vapor ) vapor pressure

∂Mm m +∇·b Fm ) qm + n˙Mweight ∂t

(32)

Mm ) φSAFAXAm + φSGFGXGm

(33)

b Fm ) XAmb FA + XGmb FG + b J Gm

(34)

where

b J Gm ) - b J Gw ) φSG(φ1/3SG7/3)DGwFG∇XGw ) φSG(τG)DGwFG∇XGw (35) qm ) XAmqA + XGmqG

(36)

The continuity equation for the water component, k ) w, is ∂Mw w +∇·b Fw ) qw + n˙NHMweight ∂t

(37)

Mw ) φSAFAXAw + φSGFGXGw + φSIFIXIw

(38)

b Fw ) XAwb FA + XGwb FG + b J Gw

(39)

qw ) XAwqA + XGwqG

(40)

where

Greek Letters ∆E ) activation energy (J mol-1 or J kg-1) ∆HH0 ) latent heat of hydrate formation or dissociation (J mol-1 or J kg-1) ∆HI0 ) latent heat of ice formation or dissociation (J mol-1 or J kg-1) ∆t ) time step (s) µA ) aqueous-phase viscosity (Pa s) µG ) gas-phase viscosity (Pa s) Fβ ) density of phase β (kg m-3) σ0 ) Stefan-Boltzmann constant (5.6687 × 10-8 J m-2 K-4) φ ) porosity of the porous media

5241

The continuity equation for the hydrate component, k ) h, is ∂Mh h ) -n˙Mweight ∂t

(41)

The energy balance equation, for k ) θ, is ∂Mθ +∇·b F θ ) qθ ∂t

(42)

5242

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

∑

Mθ ) (1 - φ)FRHR +

φSβFβHβ

(43)

β)A,G,H,I

b Fθ ) -Kave∇T + fσσ0∇T4 +

∑

Hβb Fβ

(44)

β)A,G

qθ )

∑

Hβqβ + n˙∆HH0

(45)

β)A,G

The new governing eqs 32, 37, 41, and 42 are discretized in time using an implicit method to maintain the numerical stability in the simulation as follows ∂Mkn 1 ) ∂t Vn ∆t

k Mn,t+∆t

Mkn,t

∑

1 ) Vn

k AnmFnm + qkn

(46)

f1(x1, x2, ..., xn) ) 0 f2(x1, x2, ..., xn) ) 0 l fn(x1, x2, ..., xn) ) 0

f(x + ∆x) ) f(x) + J(x)∆x

∑

+

k qn,t+∆t

∑

k k AnmFnm,t+∆t - ∆tqn,t+∆t )0

nm

(48) For the kinetic model, eq 48 actually represents the four equations, for the components methane, water, hydrate, and energy. The next step is to solve eq 48 to obtain the primary variables at a new time, t + ∆t. This is discussed in the next section. Newton-Raphson Method The Newton-Raphson algorithm is one of the best-known methods of finding roots because of its simplicity and speed.38 The Newton-Raphson formula can be derived from the Taylor series expansion of function f(x) around x f(xi+1) ) f(xi) + f'(xi)(xi+1 - xi) + O(xi+1 - xi)2

(49)

If xi+1 is a root of f(x) ) 0, then f(xi+1) ) 0, and eq 49 becomes 0 ) f(xi) + f'(xi)(xi+1 - xi) + O(xi+1 - xi)2

(50)

Dropping the higher-order term in eq 50 gives the NewtonRaphson formula xi+1 ) xi -

f(xi) f'(xi)

(51)

The algorithm for the Newton-Raphson method is to repeat eq 51 starting with an initial value of x0, until the convergence criterion is reached, namely |xi+1 - xi | < ε

∂fi ∂xj

Because analytical derivation of each ∂fi/∂xj expression can be difficult or impractical, it is preferable to let the computer calculate the partial derivatives from the finite-difference approximation ∂fi fi(xj + 10-8xj) - fi(xj) ≈ ∂xj 10-8xj

∆x ) -

f(x) J(x)

(58)

The following steps constitute the Newton-Raphson method for simultaneous nonlinear equations: (1) Estimate the solution vector x. (2) Evaluate f(x). (3) Compute the Jacobian matrix J(x) from eq 57. (4) Set up the simultaneous equations in eq 58 and solve for ∆x. (5) Let xnew r xold + ∆x and repeat steps 2-5. The above process is continued until |∆x| < ε, where ε is the error tolerance and was set to 10-5 in the simulations reported here. To show the details of the Newton-Raphson method for the kinetic model, assume that, at time t, there are two elements (i and j) in the simulation as shown in Figure 8. Element i has gas and hydrate phases, and its primary variables are xi,1 (gas pressure PG), xi,2 (gas saturation SG), xi,3 (ice saturation SI), and xi,4 (temperature T). Element j has aqueous and hydrate phases, and its primary variables are xj,1 (water pressure PA), xj,2 (aqueous-phase saturation SA), xj,3 (mass fraction of methane in the aqueous phase Xm A ), and xj,4 (temperature T). Note that the primary variables might

(52)

(53)

(57)

Assuming that x is the current approximation of the solution of f(x) ) 0, and letting x + ∆x be the improved solution, the correction ∆x is found by setting f(x + ∆x) ) 0 in eq 55. The result is a set of linear equations for ∆x

where ε is the error tolerance. The default is ε ) 10-5 in HydrateResSim. Consider the n-dimensional version of the same problem, namely f(x) ) 0

(56)

(47)

nm

∆t Vn

(55)

where J(x) is the Jacobian matrix, of size n × n, made up of the partial derivatives Jij )

k AnmFnm,t+∆t

(54)

To derive the Newton-Raphson method for a system of equations, we drop the high-order terms and start with the Taylor series expansion about the point x

m

The implicit-time method means that the flux and source terms on the right-hand side of eq 47 are evaluated at time t + ∆t. We define the residual, as shown in eq 48, as the difference between the left- and right-hand sides of eq 47 k k Rn,t+∆t ) Mn,t+∆t - Mkn,t -

or, using scalar notation

Figure 8. Schematic for elements i and j.

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

not be the same at all nodes in the domain. They can be substituted during the iteration procedure within the calculation of a time step if there is a phase change.

m m Ri,t+∆t ) fi,1 ) Mi,t+∆t - Mm i,t -

∆t Vi

∑A

m imFim,t+∆t

∂fi,1 ∂xi,1 ∂fi,2 ··· ∂xi,1 ∂fi,3 ··· ∂xi,1 ∂fi,4 ∂xi,1 J ) l ∂fj,1 ∂xi,1 ∂fj,2 ··· ∂xi,1 ∂fj,3 ∂xi,1 ∂fj,4 ··· ∂xi,1 l

-

im

m m ∆tqi,t+∆t - ∆tn˙i,t+∆tMweight ) 0(59)

w w Ri,t+∆t ) fi,2 ) Mi,t+∆t - Mwi,t -

∆t Vi

∑A

w imFim,t+∆t

-

im

w w ∆tqi,t+∆t - ∆tn˙i,t+∆tNHMweight ) 0(60)

∆t Vi

∑A

θ imFim,t+∆t

-

im

θ ∆tqi,t+∆t - ∆tn˙i,t+∆t∆HH0 ) 0(62)

Similarly, we have the following equations for element j m m Rj,t+∆t ) fj,1 ) Mj,t+∆t - Mm j,t -

∆t Vj

∑A

m jmFjm,t+∆t

-

jm

m m ∆tqj,t+∆t - ∆tn˙j,t+∆tMweight ) 0(63)

w w Rj,t+∆t ) fj,2 ) Mj,t+∆t - Mwj,t -

∆t Vj

∑A

w jmFjm,t+∆t

-

jm

w w ∆tqj,t+∆t - ∆tn˙j,t+∆tNHMweight ) 0(64)

h Rj,t+∆t

) fj,3 )

h Mj,t+∆t

-

Mhj,t

θ θ Rj,t+∆t ) fj,4 ) Mj,t+∆t - Mθj,t -

+

∆t Vj

h ∆tn˙j,t+∆tMweight

∑A

θ jmFjm,t+∆t

∂fi,1 ∂xi,2 ∂fi,2 ∂xi,2 ∂fi,3 ∂xi,2 ∂fi,4 ∂xi,2

∂fi,1 ∂xi,3 ∂fi,2 ∂xi,3 ∂fi,3 ∂xi,3 ∂fi,4 ∂xi,3

∂fi,1 ∂fi,1 ··· ∂xi,4 ∂xj,1 ∂fi,2 ∂fi,2 ··· ∂xi,4 ∂xj,1 ∂fi,3 ∂fi,3 ··· ∂xi,4 ∂xj,1 ∂fi,4 ∂fi,4 ∂xi,4 ∂xj,1

∂fi,1 ∂xj,2 ∂fi,2 ∂xj,2 ∂fi,3 ∂xj,2 ∂fi,4 ∂xj,2

∂fi,1 ∂xj,3 ∂fi,2 ∂xj,3 ∂fi,3 ∂xj,3 ∂fi,4 ∂xj,3

∂fj,1 ∂xi,2 ∂fj,2 ∂xi,2 ∂fj,3 ∂xi,2 ∂fj,4 ∂xi,2

∂fj,1 ∂xi,3 ∂fj,2 ∂xi,3 ∂fj,3 ∂xi,3 ∂fj,4 ∂xi,3

∂fj,1 ∂xi,4 ∂fj,2 ∂xi,4 ∂fj,3 ∂xi,4 ∂fj,4 ∂xi,4

∂fj,1 ∂xj,1 ∂fj,2 ··· ∂xj,1 ∂fj,3 ··· ∂xj,1 ∂fj,4 ··· ∂xj,1

∂fj,1 ∂xj,2 ∂fj,2 ∂xj,2 ∂fj,3 ∂xj,2 ∂fj,4 ∂xj,2

∂fj,1 ∂fj,1 ∂xj,3 ∂xj,4 ∂fj,2 ∂fj,2 ··· ∂xj,3 ∂xj,4 ∂fj,3 ∂ j,3 ··· ∂xj,3 ∂xj,4 n> ∂fj,4 ∂fj,4 ··· ∂xj,3 ∂xj,4

···

∂fi,1 ∂xj,4 ∂fi,2 ∂xj,4 ∂fi,3 ∂xj,4 ∂fi,4 ∂xj,4

··· ···

where xi,1 is the first primary variable in element i, xi,2 is the second primary variable in element i, xi,3 is the third primary variable in element i, xi,4 is the fourth primary variable in element i, xj,1 is the first primary variable in element j, xj,2 is the second primary variable in element j, xj,3 is the third primary variable in element j, and xj,4 is the fourth primary variable in element j. fi,1 is a function not only of the primary variables its own element xi,1, xi,2, and xi,3 but also of the other element’s primary variables because of the interaction between different elements such as the flux of mass and energy transfer. The effects of the different primary variables for a given phase combination are found through the evaluation of the Jacobian matrix. For example, we use eq 57 to evaluate the derivatives of functions

h h h Ri,t+∆t ) fi,3 ) Mi,t+∆t - Mhi,t + ∆tn˙i,t+∆tMweight ) 0 (61)

θ θ Ri,t+∆t ) fi,4 ) Mi,t+∆t - Mθi,t -

(67)

l

The nonlinearities in the equations are handled with the Newton-Raphson method. We can write the system of four equations given by eq 48 simplified in the following functional form

5243

) 0 (65)

∂fi,2 fi,2(..., xi,1 + 10-8xi,1, ...) - fi,2(..., xi,1, ...) ≈ ∂xi,1 10-8xi,1

(68)

∂fi,2 fi,2(..., xi,2 + 10-8xi,2, ...) - fi,2(..., xi,2, ...) ≈ ∂xi,2 10-8xi,2

(69)

∂fi,2 fi,2(..., xi,3 + 10-8xi,3, ...) - fi,2(..., xi,3, ...) ≈ ∂xi,3 10-8xi,3

(70)

∂fi,2 fi,2(..., xi,4 + 10-8xi,4, ...) - fi,2(..., xi,4, ...) ≈ ∂xi,4 10-8xi,4

(71)

When we choose different primary variables to be solved, the partial derivative ∂fi,2/∂xi,1, ∂fi,2/∂xi,2, ∂fi,2/∂xi,3, and ∂fi,2/∂xi,4 will be different. These differences in the partial derivatives affect the solution of the governing equations.

-

jm

θ ∆tqj,t+∆t - ∆tn˙j,t+∆t∆HH0 ) 0(66)

Primary Switch Method for the Equilibrium Model In the equilibrium model, the hydrate is treated as a phase with a fixed composition, that is, the water and methane

The Jacobian matrix can be written as

Table 3. Primary Variables (PVs) in the Equilibrium Model

PV1 PV2 PV3

G

A

A+G

I+G

A+H

I+H

A+H+G

A+I+G

A+I+H

I+H+G

I+H+A+G

PG YGm T

PA XAm T

PG SA T

PG SI T

PA SA T

PA SI T

SG SA T

PG SA SG

PA SA SI

SG SI T

SG SA SI

5244

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010

molecular and weight ratios are constant in the hydrate phase. In the equilibrium model, it is assumed that the hydrate phase (h) is composed of the components methane (m) and water (w). The other assumptions remain the same as in the equilibrium model. The continuity equation for the methane component, k ) m, is

w w w Rj,t+∆t ) fj,2 ) Mj,t+∆t - Mj,t -

∆t Vj

∑A

w jmFjm,t+∆t

-

jm

w ∆tqj,t+∆t ) 0(81) θ θ θ Rj,t+∆t ) fj,3 ) Mj,t+∆t - Mj,t -

∆t Vj

∑A

θ jmFjm,t+∆t

-

jm

θ ∆tqj,t+∆t ) 0(82)

m

∂M +∇·b F m ) qm ∂t

(72)

Mm ) φSAFAXAm + φSGFGXGm + φSHFHXHm

(73)

The continuity equation for the water component, k ) w,

The Jacobian matrix can be written as l ···

is: ···

∂Mw +∇·b F w ) qw ∂t

(74)

···

∂fi,1 ∂xi,1 ∂fi,2 ∂xi,1 ∂fi,3 ∂xi,1

∂fi,1 ∂xi,2 ∂fi,2 ∂xi,2 ∂fi,3 ∂xi,2

∂fi,1 ∂xi,3 ∂fi,2 ∂xi,3 ∂fi,3 ∂xi,3

∂fj,1 ∂xi,1 ∂fj,2 ∂xi,1 ∂fj,3 ∂xi,1

∂fj,1 ∂xi,2 ∂fj,2 ∂xi,2 ∂fj,3 ∂xi,2

∂fj,1 ∂xi,3 ∂fj,2 ∂xi,3 ∂fj,3 ∂xi,3

··· ··· ···

∂fi,1 ∂xj,1 ∂fi,2 ∂xj,1 ∂fi,3 ∂xj,1

∂fi,1 ∂xj,2 ∂fi,2 ∂xj,2 ∂fi,3 ∂xj,2

∂fi,1 ∂xj,3 ∂fi,2 ∂xj,3 ∂fi,3 ∂xj,3

∂fj,1 ∂xj,1 ∂fj,2 ∂xj,1 ∂fj,3 ∂xj,1

∂fj,1 ∂xj,2 ∂fj,2 ∂xj,2 ∂fj,3 ∂xj,2

∂fj,1 ∂xj,3 ∂fj,2 ∂xj,3 ∂fj,3 ∂xj,3

··· ··· ···

(83)

J ) l

where Mw ) φSAFAXAw + φSGFGXGw + φSHFHXHw + φSIFIXIw (75)

···

The energy balance equation, k ) θ, is

···

∂Mθ +∇·b F θ ) qθ ∂t

(76)

The following comments pertain to discretization and the use of the Newton-Raphson method. There are three partial differential equations in the equilibrium model. Thus, three primary variables are needed in the equilibrium model. Table 3 lists the primary variables in the equilibrium model used in HydrateResSim for different combinations of phases. To show the details of the Newton-Raphson method, we assume that, at time t, there are two elements (i and j) in our simulation, as shown in Figure 8. Element i has aqueous and gas phases (A + G), shown in the fourth column in Table 3, and its primary variables are xi,1 (pressure PG), xi,2 (aqueousphase saturation SA), and xi,3 (temperature T). Element j has aqueous and hydrate phases (A + H), as shown in the sixth column in Table 3, and its primary variables are xj,1 (pressure PA), xj,2 (aqueous-phase saturation SA), and xj,3 (temperature T). Therefore m m Ri,t+∆t ) fi,1 ) Mi,t+∆t - Mm i,t -

∆t Vi

∑A

m imFim,t+∆t

-

im

m ∆tqi,t+∆t ) 0(77)

w w Ri,t+∆t ) fi,2 ) Mi,t+∆t - Mwi,t -

∆t Vi

∑A

w imFim,t+∆t

-

im

w ∆tqi,t+∆t ) 0(78)

θ θ Ri,t+∆t ) fi,3 ) Mi,t+∆t - Mθi,t -

∆t Vi

∑A

θ imFim,t+∆t

-

im

θ ∆tqi,t+∆t ) 0(79)

Similarly, we have the following equations for element j m m Rj,t+∆t ) fj,1 ) Mj,t+∆t - Mm j,t -

∆t Vj

∑A

···

m jmFjm,t+∆t

-

jm

m ∆tqj,t+∆t ) 0(80)

··· ··· ···

··· ··· ···

l

Literature Cited (1) Gidaspow, D. Multiphase Flow and Fluidization: Continuum Kinetic Theory Description; Academic Press: New York, 1994. (2) Peaceman, D. W. Fundamentals of Numerical ReserVoir Simulation; Elsevier: New York, 1977. (3) Boswell, R. Resource potential of methane hydrate coming into focus. J. Pet. Sci. Eng. 2007, 56, 9–13. (4) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press: Boca Raton, FL, 2007. (5) Moridis, G. J.; Collett, T. S.; 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. SPE ReserVoir EVal. Eng. 2009, 12, 745–771. (6) Max, M. D.; Johnson, A. H.; Dillon, W. P. Economic Geology of Natural Gas Hydrate; Springer: Dordrecht, The Netherlands, 2006. (7) Goel, N. In situ methane hydrate dissociation with carbon dioxide sequestration: Current knowledge and issues. J. Pet. Sci. Eng. 2006, 51, 169–184. (8) Holder, G. D.; Angert, P. F. Simulation of gas production from a reservoir containing both gas hydrates and free natural gas. Presented at the 1982 Annual Fall Technical Conference and Exhibition, New Orleans, LA, Sep 26-29, 1982; Paper SPE 11105. (9) Selim, M. S.; Sloan, E. D. Heat and mass transfer during the dissociation of hydrate in porous media. AIChE J. 1989, 35.6, 1049–1052. (10) Selim, M. S.; Sloan, E. D. Hydrate Dissociation in Sediment. SPE ReserVoir Eng. 1990, 5, 245–251. (11) Yousif, M. H.; Abass, H. H.; Selim, M. S.; Sloan, E. D. Experimental and theoretical investigation of methane gas hydrate dissociation in porous media. SPE ReserVoir Eng. 1991, 69–76. (12) Islam, M. R. A new recovery technique for gas production from Alaskan gas hydrates. J. Pet. Sci. Eng. 1994, 11, 267–281. (13) Makogon, Y. F. Hydrates of Hydrocarbons; PennWell Publishing Company: Tulsa, OK, 1997. (14) Moridis, G. J.; Apps, J.; Pruess, K.; Myer, L. EOSHYDR: A TOUGH2 Module for CH4-Hydrate Release and Flow in the Subsurface; Report LBNL-42386; Lawrence Berkeley National Laboratory: Berkeley, CA, 1998. (15) Moridis, G. J. Numerical Simulation Studies of Thermally-Induced Gas Production From Hydrate Accumulations With No Free Gas Zones at the Mallik Site, Mackenzie Delta, Canada. Presented at the SPE Asia Pacific Oil and Gas Conference and Exhibition, 8-10 October 2002; Paper SPE 77861.

Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 (16) Moridis, G. J. Numerical studies of gas production from class 2 and class 3 hydrate accumulations at the mallik site, Mackenzie Delta, Canada. SPE ReserVoir EVal. Eng. 2004, 7, 175–183. (17) Moridis, G. J.; Kowalsky, M. B.; Pruess, K. HydrateResSim User’s Manual: A Numerical Simulator for Modeling the BehaVior of Hydrates in Geologic Media; www.netl.doe.gov. (18) Moridis, G. J.; Sloan, E. D. Gas production potential of disperse low-saturation hydrate accumulations in oceanic sediments. Energy ConVers. Manage. 2007, 48, 1834–1849. (19) Moridis, G. J.; Kowalsky, M. B.; Pruess, K. TOUGH+HYDRATE User’s Manual: A Code for the Simulation of System BehaVior in HydrateBearing Geologic Media; www.esd.lbl.gov. (20) Ahmadi, G.; Ji, C.; Smith, D. H. Numerical solution for natural gas production from methane hydrate dissociation. J. Pet. Sci. Eng. 2004, 41, 269–285. (21) Ahmadi, G.; Ji, C.; Smith, D. H. Production of natural gas from methane hydrate by a constant downhole pressure well. Energy ConVers. Manage. 2007, 48, 2053–2068. (22) Ahmadi, G.; Ji, C.; Smith, D. H. Natural gas production from hydrate dissociation: An axisymmetric model. J. Pet. Sci. Eng. 2007, 58, 245–258. (23) Clarke, M. A.; Bishnoi, P. R. Measuring and modeling the rate of decomposition of gas hydrates formed from mixtures of methane and ethane. Chem. Eng. Sci. 2001, 56, 4715–4724. (24) Masuda, Y.; Fujinaga, Y.; Naganawa, S.; Fujita, K.; Sato, K.; Hayashi, Y. Modeling and Experimental Studies on Dissociation of Methane Gas Hydrates in Berea Sandstone Cores. In Proceedings of the Third International Conference on Gas Hydrates; 1999. Salt Lake City, Utah, USA. (25) Masuda, Y.; Kurihara, M.; Ohuchi H.; Sato T. A Field-Scale Simulation Study on Gas Productivity of Formations Containing Gas Hydrates. In Proceedings of the Fourth International Conference on Gas Hydrates; 2002. Yokohama, Japan, pp 40–46. (26) Ji, C.; Ahmadi, G.; Smith, D. H. Natural gas production from hydrate decomposition by depressurization. Chem. Eng. Sci. 2001, 56, 5801– 5814. (27) Ji, C.; Ahmadi, G.; Smith, D. H. Constant rate natural gas production from a well in a hydrate reservoir. Energy ConVers. Manage. 2003, 44, 2403–2423.

5245

(28) Sun, X.; Nanchary, N.; Mohanty, K. K. 1-D modeling of hydrate depressurization in porous media. Transp. Porous Media 2005, 58, 315– 338. (29) Sun, X.; Mohanty, K. K. Kinetic simulation of methane hydrate formation and dissociation in porous media. Chem. Eng. Sci. 2006, 61, 3476– 3495. (30) Tsimpanogiannis, I. N.; Lichtner, P. C. Parametric study of methane hydrate dissociation in oceanic sediments driven by thermal stimulation. J. Pet. Sci. Eng. 2007, 56, 165–175. (31) Castaldi, M. J.; Zhou, P. C.; Yegulalp, T. M. Down-hole combustion method for gas production from methane hydrates. J. Pet. Sci. Eng. 2007, 56, 176–185. (32) Kowalsky, M. B.; Moridis, G. J. Comparison of kinetic and equilibrium reaction models in simulating gas hydrate behavior in porous media. Energy ConVers. Manage. 2007, 48, 1850–1863. (33) Liu, Y.; Strumendo, M.; Arastoopour, H. Numerical Simulation of Methane Production from a Methane Hydrate Formation. Ind. Eng. Chem. Res. 2008, 47, 2817–2828. (34) Liu, Y.; Strumendo, M.; Arastoopour, H. Simulation of Methane Production from Hydrates by Depressurization and Thermal Stimulation. Ind. Eng. Chem. Res. 2009, 48, 2451–2464. (35) Nazridoust, K.; Ahmadi, G. Computational modeling of methane hydrate dissociation in a sandstone core. Chem. Eng. Sci. 2007, 62, 6155– 6177. (36) Kim, H. C.; Bishnoi, P. R.; Heidemann, R. A.; Rizvi, S. S. H. Kinetics of methane hydrate decomposition. Chem. Eng. Sci. 1987, 42, 1645–1653. (37) Gamwo, I. K.; Myshakin, E. M.; Zhang, W.; Warzinski, R. P. Simulation of Methane Production in a Laboratory-Scale Reactor Containing Hydrate-Bearing Porous Medium. Prepr. Pap.-Am. Chem. Soc., DiV. Fuel Chem. 2008, 53 (1), 454–455. (38) Kiusalaas, J. Numerical Methods in Engineering with Matlab; Cambridge University Press: New York, 2005.

ReceiVed for reView September 15, 2009 ReVised manuscript receiVed December 30, 2009 Accepted January 5, 2010 IE901452V