Ind. Eng. Chem. Res. 2009, 48, 2451–2464
2451
Simulation of Methane Production from Hydrates by Depressurization and Thermal Stimulation Yong Liu, Matteo Strumendo, and Hamid Arastoopour* Department of Chemical and Biological Engineering, Illinois Institute of Technology, 10 W. 33rd Street, Chicago, Illinois 60616
Recently methane hydrates have attracted attention due to their large quantity on the earth and their potential as a new resource of energy. This paper describes a one-dimensional mathematical model and numerical simulation of methane hydrate dissociation in hydrate reserves by both depressurization and thermal stimulation using a onedimensional radial flow system (axisymmetric reservoir). A moving front that separates the hydrate reserve into two zones is included in this model. A numerical coordinate transformation method was used to solve the moving boundary problem. The partial differential equations were discretized into ordinary differential equations using the method of lines. Our simulations showed that the moving front location and the gas flow rate production are strong functions of the well pressure and reservoir temperature. The impermeable boundary condition at the reservoir results in very low temperature at the moving front and the formation of ice. The formation of ice, which plugs the pore volume for the gas to flow, should be avoided. Compared with a stationary water phase model, our simulations showed that the assumption of a stationary water phase overpredicts the location of the moving front and the dissociation temperature at the moving front and underpredicts the gas flow rate. The thermal stimulation using constant temperature at the well method using a single well was found to have a limited effect on gas production compared to gas production due to depressurization. 1. Introduction Natural gas hydrates are solid molecular compounds composed of water with natural gas; there are huge amounts of methane gas trapped in hydrate reservoirs (Makogon,1 Ahmadi et al.,2 Sloan and Koh3). Three main methods exist to dissociate the hydrate: depressurization, inhibitor stimulation, and thermal stimulation (Ji et al.4). Mathematical modeling and simulation of the hydrate dissociation process in the hydrate deposits is essential in order to evaluate the rates of natural gas production attainable. Selim and Sloan5 developed a one-dimensional model considering the convective and conductive heat transfer in the gas hydrate reserve. They obtained an analytical solution under the assumption that a moving front separates the reserve into two zones: the dissociated gas zone and the hydrate zone. The water phase is stationary after hydrate dissociation in porous media. Ahmadi et al.2 developed a numerical method to solve a onedimensional hydrate dissociation model with both conductive and convective heat transfer and energy balance at the moving front. Ahmadi et al.6,7 extended the mathematical model into the cylindrical coordinate system. Nazridoust and Ahmadi8 developed a user-defined function that can be used by the commercial software FLUENT and solved the methane hydrate dissociation problem, and their results compared well with the experimental data by Masuda et al.9 Besides the mathematical modeling based on the moving boundary model that assumes that the hydrate dissociates only in the moving front, there is another category that is based on the assumption that the hydrate dissociates throughout the entire reservoir. On the basis of the general-purpose multiphase, multicomponent simulator TOUGH2 of fluid flow and heat transport in geologic media, Moridis et al.,10 Alp et al.,11 and Moridis and Sloan12 simulated the hydrate dissociation using a reaction model. * To whom correspondence should be addressed. E-mail: arastoopour@ iit.edu.
Liu et al.13 developed a one-dimensional two-phase flow model and numerical solutions coupled with the methane hydrate dissociation moving front concept using a Cartesian coordinate system. In our earlier studies, the mathematical model used by Ahmadi et al.7 was modified to account for the water movement (Liu et al.13). Here, the numerical techniques of Liu et al.13 are used to simulate the moving boundary problem in the axisymmetric system. The effects of the assumption of a stationary water phase on the moving front speed and the gas flow rate are discussed. In addition the effect of thermal stimulation on the methane gas production, the effect of the effective gas permeability in the undissociated zone, and the effect of the boundary condition at the reservoir are presented. 2. Hydrate Dissociation Model We simulated the case in which a single vertical production well of radius r0 ) 0.1 m is drilled into a hydrate reservoir. In the reservoir the hydrate coexists with methane gas. The depressurization method and combined depressurization and thermal stimulation methods were used to dissociate the methane hydrate. The hydrate decomposition model may be expressed as (Max et al.14) 164 m3 (CH4)STP + 0.8 m3 (H2O)water S 1.0 m3 (hydrate) (1) The hydrate decomposes into methane gas and water when the pressure decreases or the temperature increases. Figure 1 shows that a reservoir in which hydrate coexists with gas. Initially, the gas and the hydrate coexist at a pressure of Pe (reservoir pressure) and at a temperature of Te (reservoir temperature). Pe is higher than the thermodynamic equilibrium pressure of the hydrate at Te. The hydrate near the well becomes unstable and begins to dissociate into water and gas near the well, as soon as a single vertical production well is drilled into
10.1021/ie8005275 CCC: $40.75 2009 American Chemical Society Published on Web 10/29/2008
2452 Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009
Figure 1. Schematic for hydrate dissociation in an axisymmetric system.
the reservoir and the well pressure is decreased to a value less than the thermodynamic equilibrium pressure at Te or the temperature at the well is increased due to the thermal stimulation. The hydrates dissociate in a small region, which is viewed as a moving front. This moving front separates the entire reservoir into two regions: the dissociated region (zone 1 in Figure 1) consisting of sand, gas, and water, and the undissociated region (zone 2 in Figure 1) consisting of sand, gas, and hydrate. The moving front expands with time. The natural gas moves toward the production well due to the reverse pressure drop toward the production well. 3. Mathematical Model In the moving front model for hydrate dissociation, the hydrates are assumed to dissociate in a small region and this small region is assumed to be a moving front. The moving front separates the entire simulation domain into two zones: zone 1 is the dissociated zone consisting of sand, gas, and water; zone 2 is the undissociated zone consisting of sand, gas, and hydrate as mentioned. The governing equations for hydrate dissociation in porous media include the following equations: the conservation of mass, momentum, and energy in both zones and at the moving front. The governing equations used by Liu et al.13 were transformed into the a radial coordinate system. A total length of R0 ) 100 m of hydrate reservoir was simulated in this study. Continuity equation for water in zone 1: ∂(φSwFw) 1 ∂(rFwUw) + )0 ∂t r ∂r
ΚabsΚrw ∂Pw µw ∂r
Κrw ) Sw4
∂(φSg1Fg1) 1 ∂(rFg1Ug1) + )0 ∂t r ∂r
(4)
Darcy’s law for the superficial gas velocity in zone 1: Ug1 ) -
ΚabsΚrg ∂Pg1 µg ∂r
Κrg ) (1 - Sw)2(1 - Sw2)
(5)
Capillary pressure: Pc ) Pg1 - Pw ≈ 0
(6)
The capillary pressure is the pressure difference between the gas pressure and the water pressure. The capillary pressure is caused by the surface tension and the capillary force between the gas and the water phase. For methane hydrate dissociation in porous media, the capillary pressure is very small when compared with the gas pressure (Yousif et al.16); thus, the capillary pressure effect is neglected in the simulations. The water and gas saturations have the following relationship in zone 1: Sw + Sg1 ) 1
(7)
The density of gas is a function of temperature and pressure: Fg )
PgM ZRT
(8)
Continuity equation for the gas phase in zone 2: (2)
Darcy’s law for the superficial water velocity in zone 1 (Corey15): Uw ) -
Continuity equation for gas in zone 1:
(3)
∂(φSg2Fg2) 1 ∂(rFg2Ug2) + )0 ∂t r ∂r
(9)
Darcy’s law for the superficial gas velocity in zone 2: Ug2 ) -
Κeff2 ∂Pg2 µg ∂r
(10)
Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009 2453
The formula to calculate effective permeability is the same as the one used in Liu et al.13 Energy conservation equation in zone 1 (Masuda et al.9): ∂(φSwFwCpwT1 + φSg1Fg1CpgT1 + (1 - φ)FsCpsT1) + ∂t 1 ∂(rFg1CpgUg1T1 + rFwCpwUwT1) ) r ∂r 1 ∂(rkgφSg ∇T1 + rkwφSw ∇T1 + rks(1 - φ) ∇T1) (11) r ∂r Energy conservation equation in zone 2: ∂(φShFhCphT2 + φSg2Fg2CpgT2 + (1 - φ)FsCpsT2) + ∂t 1 ∂(rFg2CpgUg2T2) ) r ∂r 1 ∂(rkgφSg ∇T2 + rkhφSh ∇T2 + rks(1 - φ) ∇T2) (12) r ∂r Mass balance for gas at the moving front:
T1(r0, t) ) Twell for thermal stimulation combined with depressurization (22) S(0) ) r0
(23)
Pwell and Twell are the gas pressure and temperature at the production well. R0 is the total length of simulation. Pe is the gas pressure of the reservoir which is a constant and equals to 15 MPa. Te is the temperature of the reservoir and equals to 287 K. PD and TD are the dissociation pressure and dissociation temperature at the moving front, respectively. Both are a function of time and are calculated numerically during the simulations. The location of the moving front S(t) is a function of time and its initial location is at the production well. The conservation of water at the hydrate dissociation moving front results in the boundary condition for the water saturation. Sw )
(1 - ε)FhSh Fw
(24)
In equation 24, Sh is the hydrate saturation in zone 2 and remains constant during the simulations, ε is the mass fraction of methane per unit mass of hydrate and equals to 0.129, Fh is the density of hydrate, and Fw is the density of water. 4. Coordinate Transformation Method
dS 2πrφ (ShεFh + Sg2Fg2 - Sg1Fg1) ) 2πrFg2Ug2 - 2πrFg1Ug1 dt (13) Energy balance at the moving front: dS ) 2πrFg1CpgT1Ug1 + 2πrFwCpwT1Uw dt 2πr(kg1φSg1 + kwφSw + ks(1 - φ)) ∇T1 + 2πr(kg2φSg2 + khφSh + ks(1 - φ)) ∇T2 - 2πrFg2CpgT2Ug2 (14)
2πr∆HFhφSh
The main idea is to transform the governing equations into a new coordinate system so that the moving boundary’s position is fixed in the new coordinate system. Wang et al.17 applied the numerical coordinate transformation to the solidification of binary alloys. Liu et al.13 applied the coordinate transformation method to the numerical simulation of methane hydrate dissociation in a Cartesian coordinate system. The coordinate transformation method used by Liu et al.13 was modified to be used in the radial flow system, according to the following equations. zone 1
Thermodynamic equilibrium at the front:
ξ)
r - r0 S(t) - r0
0
