Numerical Simulation of Methane Production from ... - ACS Publications

---

Ind. Eng. Chem. Res. 2008, 47, 2817-2828

2817

GENERAL RESEARCH Numerical Simulation of Methane Production from a Methane Hydrate Formation Yong Liu, Matteo Strumendo, and Hamid Arastoopour* Department of Chemical and Biological Engineering, Illinois Institute of Technology, 10 West 33rd Street, Chicago, Illinois 60616

This paper describes a one-dimensional model for hydrate dissociation in porous media by the depressurization method. A moving boundary, which separates the total simulation zone into two zones, is used. The governing equations consider the convective-conductive heat transfer and mass transfer in the gas and hydrate zones together with the energy balance at the moving front. These equations were transformed into a new coordinate system using a coordinate transformation method. The numerical method of lines was used to discretize the governing equations after coordinate transformation. Distributions of temperature and pressure for different well pressure and reservoir temperature are presented. The speed of the moving front and the gas production rate were shown to be strong functions of the well pressure and the absolute permeability of the porous media. Our simulations also showed that the assumption of stationary water phase, underpredicts gas production and overpredicts the speed of the moving front. 1. Introduction Natural gas hydrates are solid molecular compounds composed of water and a large amount of methane gas is trapped in the reservoirs.1-3 Methane is a less carbon-intensive fuel than other hydrocarbons, such as oil or coal. The combustion of methane yields 44% less CO2 than coal, per unit energy release, and 29% less CO2 than oil.2 Natural gas hydrates occur in two zones: in permafrost and under the sea floor. Because the hydrate density is smaller than the density of seawater, the hydrate is cemented with the sediment of the sea floor to be stable.2 There are three main methods to dissociate the hydrate: depressurization,4-9 inhibitor stimulation, and thermal stimulation.10-13 Mathematical modeling and simulation of the hydrate dissociation process and of the transport phenomena associated with the flow of natural gas in the hydrate deposits is essential in order to evaluate the rate of natural gas production attainable. Consequently, it is also necessary to determine the importance of the hydrate reserves in the panel of the future energy supplies. Selim and Sloan14 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. They also assumed that no gas phase exists in the hydrate zone and the water phase is stationary after hydrate dissociation in porous media. Makogon15 proposed a model with an analytical solution considering the adiabatic and throttling effects of methane gas, together with the gas convective energy, with the assumption that both gas phase and hydrate phase exist in the hydrate zone. Ji et al.7 calculated the undetermined values in Makogon’s analytical solution and applied it to given sets of reservoirs with different operational parameters. Both Makogon and Ji’s models neglected the energy balance at the moving front.

Ahmadi et al.4 solved a one-dimensional model numerically with both conductive and convective heat transfers and an energy balance at the moving front with the assumption of stationary water phases in the dissociated zone. They also considered only the convective heat transfer in the dissociated gas zone in the energy balance equation. More recently, Nazridoust and Ahmadi12 simulated the hydrate dissociation using commercial software FLUENT. Our model followed the approach of Ahmadi.4 We modified the governing equations to account for the water movement. The energy balance equations include heat conduction in sand that composed the porous media; and the convective heat transfer due to the movement of the gas and water after hydrate dissociation. A coordinate transformation method was applied, to avoid the calculation of the moving front explicitly. 2. Hydrate Dissociation Model Depending on the physical property of the hydrate reservoir, the gas hydrate may exist by itself or coexist with the pressurized gas. We simulated the case that a well that is drilled into a hydrate reservoir coexists with methane gas. The depressurization method was used to dissociate the methane hydrate. The depressurization method is the least expensive method of hydrate dissociation and most feasible when the gas hydrate-bearing sands have a large free gas.2 The hydrate decomposition model may be expressed as2

164 m3(CH4)STP + 0.8 m3(H2O)water T 1.0 m3(Hydrate) (1) The hydrate decomposes into methane gas and water when the pressure decreases or the temperature rises. Figure 1 shows a schematic of a reservoir of hydrate coexisting with gas. Initially, the gas and hydrate coexist at pressure of Pe (reservoir pressure) and at temperature of Te (reservoir temperature). Pe is higher than the thermodynamic

10.1021/ie071398b CCC: $40.75 © 2008 American Chemical Society Published on Web 03/19/2008

2818

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Figure 1. Schematic for one-dimensional hydrate dissociation in porous media.

equilibrium pressure of 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 production well is drilled into the reservoir and the well pressure is decreased to a value less than the thermodynamic equilibrium pressure for hydrate at the reservoir temperature Te. 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) consists of sand, gas, and water, while the undissociated region (zone 2 in Figure 1) consists of sand, gas, and hydrate. The moving front expands as time passed. The natural gas moves toward the production well due to the reverse pressure drop toward the production well. 3. Governing Equations The heat and fluid flow, through a porous medium saturated with fluid (liquid or gas), is a function of several parameters such as the following: the porosity, the permeability of the porous medium, and the volume-averaged properties of the fluid flowing through the porous medium.16 In the derivation of the governing equations for hydrate dissociation in porous media, the porosity of the sand is assumed to be constant and the sand is stationary. The modeling of hydrate dissociation in porous media includes the conservation of mass, momentum and energy in both zones and at the moving front. We used the following governing equations to simulate the methane hydrate dissociation in porous media. These governing equations include modified equations of Ahmadi et al.4 for a two phase flow in porous media. We used the relative permeability concept to account for the two phase flow in zone 1.17-19 A total length of L ) 100 m of hydrate reservoir is simulated in this study. Continuity Equation for the Water Phase in Zone 1.

∂(φSwFw) ∂(FwUw) + )0 ∂t ∂x

(2)

Here φ is the porosity of the porous media and is assumed to be constant during hydrate dissociation. Fw is the constant water density.

Darcy’s Law for Water Velocity in Zone 1.

Uw ) -

ΚabsΚrw ∂Pw µw ∂x

Κrw ) Sw4

(3)

Here Κabs is the absolute permeability of the porous media which is a constant and equals to 4.2 millidarcy. Κrw is the relative permeability of water which is a function of the water saturation.8 Sw is water saturation and is a variable. µw is the viscosity of water and is assumed to be constant. Continuity Equation for the Gas Phase in Zone 1.

∂(φSg1Fg1) ∂(Fg1Ug1) + )0 ∂t ∂x

(4)

Darcy’s Law for Gas Velocity in Zone 1.

Ug1 ) -

ΚabsΚrg ∂Pg1 µg ∂x

Κrg ) (1 - Sw)2(1 - Sw2)

(5)

Here µg is the viscosity of gas and is assumed to be constant during the hydrate dissociation. The absolute permeability Κabs is a constant property of the porous medium. Κrg is the relative permeability of gas in zone 1 which is a function of the water saturation.19 In eqs 3 and 5, the gas pressure and the water pressure is related by

Pc ) Pg1 - Pw

(6)

In our simulations, the capillary pressure is assumed to be negligible because the capillary pressure is very small compared to the reservoir pressure.20 The relationship between the water saturation and the gas saturation may be written as

Sw + Sg1 ) 1

(7)

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2819

The density of gas is a function of temperature and pressure:

PgM ZRT

Fg )

(8)

Here M is the molecular weight of methane and R is the universal gas constant. The compressibility Z was calculated using the Redlich/Kwong equation of state.21 The calculated value of Z is 0.8838 at a pressure of 6 MPa and a temperature of 287 K. Although the compressibility Z is actually a function of temperature and pressure, it changes slightly within our pressure and temperature range. Thus, we assumed that the compressibility is constant during our simulations. Continuity Equation for the Gas Phase in Zone 2.

∂(φSg2Fg2) ∂(Fg2Ug2) + )0 ∂t ∂x

(9)

Darcy’s Law for the Superficial Gas Velocity in Zone 2

Ug2 ) -

Κeff2 ∂Pg2 µg ∂x

(10)

The Κeff2 is the effective permeability of gas with the presence of hydrate in the porous media. In zone 2, the effective permeability Κeff2 for gas is smaller than the absolute permeability because some pore space for gas flow is occupied by the stationary hydrate phase. We used the following relationship between the absolute permeability and the effective permeability by Kozeny-Carman.18

Κeff ) Κabs

the hydrate phase is stationary. Sh is the constant hydrate saturation in zone 2. Mass Balance for Gas at the Moving Front.

dS φ (βFh + (1 - β)FgD - (1 - RD)FgD) ) dt Fg2Ug2 - Fg1Ug1 (14) where  is the mass fraction of gas in methane hydrate with constant value of 0.129. The hydrate dissociation happens only at the moving front and the hydrate dissociation is an endothermic reaction. Energy Balance at the Moving Front.

∆HFhφSh

dS ) Fg1CpgT1Ug1 + FwCpwT1Uw - (kgφSg + dt

kwφSw + ks(1 - φ)) 3T1 + (kgφSg + khφSh + ks(1 - φ)) 3T2 - Fg2CpgT2Ug2 (15)

Sg3(1 - φ)2

(11)

(1 - φSg)2

Energy Conservation Equation in Zone

Figure 2. Schematic for coordinate transformation method.

1.22

∂(φSwFwCpwT1 + φSg1Fg1CpgT1 + (1 - φ)FsCpsT1) + ∂t ∂(Fg1CpUg1T1 + FwCpwUwT1) ) ∂x ∂(kgφSg∇T1 + kwφSw∇T1 + ks(1 - φ)∇T1) (12) ∂x The first term in eq 12 is the accumulation of energy for water, gas, and sand in terms of heat capacity and temperature in zone 1. The second term in eq 12 is the heat flux due to convection in the porous media which includes two terms: convective heat flux of gas and convective heat flux of water. The right-hand side of eq 12 is the conductive heat flux of gas, water, and sand. Energy Conservation Equation in Zone 2.

∂(φShFhCphT2 + φSg2Fg2CpgT2 + (1 - φ)FsCpsT2) + ∂t ∂(Fg2CpgUg2T2) ) ∂x ∂(kgφSg∇T2 + khφSh∇T2 + ks(1 - φ)∇T2) (13) ∂x The difference between eq 12 and eq 13 is that eq 13 does not have the convective heat flux term of water since we assumed

The right-hand side of eq 15 has several heat flux terms. Fg1CpgT1Ug1 + FwCpwT1Uw is the convective heat flux of gas and water from the moving front to zone 1. (kgφSg + kwφSw + ks(1 - φ))∇T1 is the conductive heat flux of gas, water and sand from the moving dissociation front to zone 1. (kgφSg + khφSh + ks(1 - φ))∇T2 is the conductive heat flux of gas, hydrate and sand from zone 2 to the moving front while Fg2CpgT2Ug2 is the convective heat flux of gas from zone 2 to the dissociation front. ∆H is the latent heat for hydrate dissociation and is assumed to be constant and equal to 458 kJ/kg.23,24 The gas and the hydrate at the moving front are assumed to be in thermodynamic equilibrium at the front and have the following relationship:4,7

log10 (PD) ) 0.0342 (TD - 273.15) + 0.0005 (TD - 273.15)2 + 6.4804 (16) The above equation is derived from the regression of the equilibrium pressure-temperature data of methane hydrate.15 Initial and Boundary Conditions. The conservation of water at the hydrate dissociation moving front results in the boundary condition for the water saturation:

FwφSw ) (1 - )FhφSh The initial and boundary conditions are listed below:1

(17)

2820

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

P1(0,t) ) Pwell

(18)

P2(x,0) ) P2(L,t) ) Pe

(19)

P1(S(t),t) ) P2(S(t),t) ) PD(TD)

(20)

T2(x,0) ) T2(L,t) ) Te

(21)

T1(S(t),t) ) T2(S(t),t) ) TD

(22)

∂T1(0,t) )0 ∂x

(23)

S(0) ) 0

(24)

where Sh is the hydrate saturation in zone 2 and remains constant during the simulation. Pwell is the gas pressure at the production well. L is the total length of simulation. Pe is the gas pressure of the reservoir. Te is the temperature of the reservoir. PD and TD are the dissociation pressure and dissociation temperature at the moving front, respectively. The location of the moving front S(t) is a function of time and its initial location is at the production well. 4. Numerical Method The position of the moving boundary is one of the unknown variables and must be calculated numerically. Ahmadi et al.4 assumed that the dissociation moving front exists initially at a 1 m long gas zone and solved the governing equations for 20 h which is estimated time for forming a 1 m long gas zone. After this initial 20 h, the location of the moving front expands with time. We avoided the limitation of assuming the formation of 1 m long gas zone by using a front fixing method for moving boundary problems. In this section, a coordinate transformation is introduced to solve the governing equations. 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.25 has applied the numerical coordinate transformation to the solidification of binary alloys. Our approach for coordinate transformation was the same as Wang. In the following, the details of the coordinate transformation are discussed: As shown in Figure 2, we assumed there is a front S(t) which separates the total domain into two regions and defined two new coordinates.

x 0 < ξ < 1 0 < x < S(t) S(t)

(25)

x - S(t) 0 < η < 1 S(t) < x < L L - S(t)

(26)

zone 1 ξ ) zone 2 η )

where L is the total length of the reservoir for simulation. Following above coordinate transformation, we transformed the governing equations into the new coordinate from the old coordinate system. In the new coordinate systems, the moving front is always at ξ ) 1 and η ) 0. The dependent variables become function of ξ and η.

U(x,t) f U(ξ,t) U(x,t) f U(η,t)

(27)

After the governing equations were transformed into the new coordinate system (see Appendix A), we then discretized these governing equations into ordinary differential equations (ODEs) using the method of lines.26 The equations were discretized in space only and kept continuous in time. We used MATLAB function ODE15s for stiff problem as the ODE solver to solve our ODEs after the discretization using the method of lines. During the discretization, a center difference for pressure term and an upwind difference for temperature techniques were used to avoid instability. The detailed governing equations after coordinate transformation are shown in Appendix A.

Figure 3. Variations of water saturation with time in the reservoir for a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy.

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2821

Figure 4. Variations of pressure distribution with time in the reservoir for a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy.

Figure 5. Variations of temperature distribution with time in the reservoir for a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy.

5. Results and Discussion 5.1. Results for Two Phase Flow Model. For all the results shown below, the total simulation length L is 100 meters. We began the dissociation by assuming that at time ) 0, the moving front already exists at the length of 10-6L, which is 0.0001 m. This is an approximation to initial condition of S(0) ) 0. Each zone was discretized into 500 intervals. The reservoir pressure was maintained at 15 MPa.

At a reservoir temperature of 287 K and well pressure of 2 MPa, Figure 3 shows the water saturation change as a function of time and space. The absolute permeability for gas and water used in this part was 4.2 millidarcy. As stated before, the moving boundary expands with time and separates the entire reservoir into two zones. Figure 3 shows that the water saturation remained constant at the moving front and then decreased toward the production well. The location of the moving front expanded

2822

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Figure 6. Comparison of water saturation distribution between different simplifications for a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy.

Figure 7. Comparison of dissociation front location between two phase flow model and stationary water phase flow model in the reservoir for a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy.

with time. After 120 days, the location of the moving front was around 16 m from the production well. The water saturation decreased with time at a higher rate at the well due to the expansion of the methane gas at the well. Ahmadi et al.1 assumed that the water saturation in zone 1 is 0.1506 and remains constant with no water flow. Figure 4 shows the variations of pressure distribution with time in the reservoir for a well pressure of 2 MPa, reservoir

temperature of 287 K, absolute permeability of 4.2 millidarcy. Figure 4 shows that the pressure decreased gradually from the reservoir to the moving front to the production well. The pressure gradient in zone 1 and zone 2 does not differ a lot since the ratio of the gas permeability in zone 1 and zone 2 was around 1.2. Figure 5 shows the variations of temperature distribution with time in the reservoir for a well pressure of 2 MPa, reservoir

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2823

Figure 8. Comparison of gas flow rate prediction using two phase flow model with the stationary water phase flow model in the reservoir for a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy.

Figure 9. Comparison of moving front location with analytical solution for a well pressure of 4 MPa, reservoir temperature of 287 K, and absolute permeability of 21 millidarcy.

temperature of 287 K, and absolute permeability of 4.2 millidarcy. Figure 5 shows that the temperature decreased slightly from the reservoir to the moving front, then decreased gradually again at a higher rate to the well. The temperature had a big sharp gradient at the moving front. This is due to the endothermic reaction of hydrate dissociation. As mentioned before, the hydrate dissociates only at the moving front. The

dissociation temperature at the moving front increased with time due to the convective and conductive heat flows from the reservoir. The temperature at the production well also increased with time due to the increase of dissociation temperature at the moving front. 5.2. Simplifying Assumptions. The governing equations may be further simplified with more assumptions such as the

2824

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Figure 10. Comparison of pressure distribution with analytical solution for a well pressure of 4 MPa, reservoir temperature of 287 K, and absolute permeability of 21 millidarcy.

Figure 11. Comparison of gas flow rate per unit area output for a reservoir temperature of 287 K and absolute permeability of 4.2 millidarcy at different well pressures.

stationary water phase. In this section, the effect of different simplifications are discussed and compared with the results of the two phase flow model. We used the Neumann boundary condition at the well for temperature in the two phase flow model as a base case. When we discussed the effect of a specific assumption, only that

specific assumption was made while all other parameters remain the same as the two phase flow model. We considered the following assumptions: (1) whether the gas density can be simplified as a function of pressure only or (2) whether the energy balance in the zone 1 is dominated by the convective heat transfer.

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2825

Figure 12. Comparison of dissociation front movement for a reservoir of 2 MPa well pressure and absolute permeability of 4.2 millidarcy at different reservoir temperatures.

Figure 6 shows the effect of the above-mentioned assumptions in the prediction of water saturation variations distribution for a well pressure of 2 MPa, reservoir temperature of 287 K and absolute permeability of 4.2 millidarcy. Figure 6 shows that neither the assumptions of gas density as a function of pressure only nor heat conduction in zone 1 affect our results. Thus, these simplifying assumptions may be made to simplify the governing equations when the depressurization method is used to dissociate methane hydrate. Figure 7 shows the comparison of the dissociation front location between the two phase flow and stationary water phase flow models in the reservoir for a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy. Figure 7 shows that the moving front differs around 1 m after 120 days. The assumption of whether the water phase is stationary is critical because the energy balance in zone 1 is mostly due to the convective heat transfer. The flow of water from the dissociation front to the well results in transfer of heat from moving front to zone 1 not only by conduction but also by convection. Thus, net heat transfer to the moving front from zone 2 is lower and in turn less heat is available for dissociation of hydrates. This lowers the speed of moving front. Therefore, the assumption of stationary water results in overpredicting the speed of the moving dissociation front. Figure 8 shows the comparison of gas flow rate as a function of time between the two phase flow and stationary water phase flow models in the reservoir with a well pressure of 2 MPa, reservoir temperature of 287 K, and absolute permeability of 4.2 millidarcy. For both models, the gas flow rate at the well increased with time due to the dissociation of hydrate. The jumps in gas flow rate in both the two phase flow model and the model with the assumption of stationary water phase show the location of the moving front. These mass flow rate jumps were due to the dissociation of hydrate at the moving front. Although Figure 7 shows that the dissociation front moved faster for the

stationary water phase flow model than the two phase model, the mass flow rate moved faster for the two phase flow model because the gas saturation increased with time in the two phase flow model while the gas saturation remained constant in the case of stationary water phase flow model. This is due to the fact that as pressure drops from moving front toward the production well gas-phase expands, but water phase remains incompressible. This results in an increase in gas saturation and in turn an increase in gas flow rate. This means that the assumption of stationary water underpredicts the gas production. 5.3. Comparison with the Analytical Solution. In this section, Makogon’s governing equations15 were solved numerically using the coordinate transformation method to validate our numerical techniques. These results were then compared with the Ahmadi et al.7 analytical solution of these equations and our two phase flow model predictions. In comparison with the analytical solution, the parameters such as the gas viscosity and absolute permeability were remained the same. Figure 9 shows the comparison of the moving front location calculated based on the two phase flow and Makogon’s equations15 for a well pressure of 4 MPa, reservoir temperature of 287 K, and absolute permeability of 21 millidarcy. The numerical solution of Makogon’s equation showed good agreement with the analytical solution by Ahmadi et al.7 Figure 9 also shows that the equations proposed by Makogon15 and later used by Ahmadi et al.7 overpredict the location of the moving front in comparison with the two phase flow model prediction. In the analytical solution, the temperature at the moving front is kept constant; thus the energy balance due to the hydrate dissociation at the moving front cannot be enforced. This explains why the speed of moving front of the two phase flow model is slower than the analytical solution.7,15 Figure 10 shows the comparison of pressure distribution predicted by our two phase flow model and numerical solution of Makogon’s equation15 using our numerical techniques with

2826

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Figure 13. Comparison of dissociation front movement for a reservoir temperature at 287 K and well pressure of 2 MPa at different absolute permeabilities.

the analytical solution. From Figure 10, the pressure gradient of the analytical solution at the well was higher than the pressure gradient of the two phase flow model. This observation enhanced our conclusion from Figure 9 that the analytical solution overpredicts the gas produced at the well. The difference between the numerical and analytical solutions of Makogon’s equation is nominal and thus validates our numerical methods. 5.4. Parametric Study. Figure 11 shows the gas output per unit area at the well at different well pressures for a reservoir temperature of 287 K and absolute permeability of 4.2 millidarcy. Figure 11 shows that a lower well pressure resulted in higher output. With the decrease of the well pressure, the dissociation front velocity increased and the gas and water flowed faster in zone 1 and gas in zone 2 also moved faster toward dissociation front. Thus, the gas in zone 2 brought more heat from the reservoir to the moving front. This is another reason for an increase in the speed of the moving front and gas output with the decrease in well pressure. Figure 12 shows the comparison of the front movement at different reservoir temperatures with the same well pressure of 2 MPa. Figure 12 shows that the higher reservoir temperature generated a higher dissociation front velocity. As mentioned before, the hydrate dissociation is an endothermic reaction. For the depressurization model, all the sensible heat for hydrate dissociation comes from the reservoir itself. Therefore, a higher reservoir temperature means that there is more sensible heat available to overcome the latent heat of hydrate dissociation and thus results higher front movement. Figure 13 shows the comparison of the moving front locations for a reservoir temperature of 287 K and well pressure of 2 MPa at different permeabilities. In Figure 13, the range of permeability is from 1.05 to 21 millidarcy. Figure 13 clearly shows that higher permeability resulted in higher front movement. The reason for this phenomenon is similar to the decrease of well pressure. Either by increasing the permeability or

decreasing the well pressure, the gas and water can move faster in zone 1 and so does the gas velocity in zone 2. This means that a higher permeability results in higher gas output. The convective heat flux also increased due to the increase of velocity of phases in both zones. Overall, our parametric studies showed that the speed of the moving front is a function of well pressure, reservoir temperature, and absolute permeability while the gas production flow rate at the well is a function of the well pressure and the absolute permeability in the range of temperatures existing in the hydrate reservoirs. 6. Conclusions One-dimensional model for hydrate dissociation in porous media was developed and solved using coordinate transformation and the numerical method of lines. The following conclusions may be derived from this study: 1. Depressurization is a feasible method for dissociation of methane hydrate. 2. The water saturation decreases gradually from the dissociation moving front to the production well in the two phase flow model. 3. Higher hydrate reservoir permeability and lower well pressure results in higher front movement and higher gas production rate at the well. 4. The assumption of stationary water phase overpredicts the location of the dissociation front and underpredicts the gas production at the well. 5. The assumption of gas density as a function of pressure can simplify the governing equation with nominal effect on the prediction of gas produced and other flow parameters. 6. In future work, we will expand our model to multidimensional domains taking into account directional permeabilities.

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2827

Acknowledgment Thanks to the Gas Technology Institute for financial support of this study and for the excellent discussions and suggestions by Mr. Iraj Salehi of the Gas Technology Institute. Appendix A. Transformed Equation under New Coordinate System

∂T2 kave2 1 - η dS ∂T2 1 × ) + 2 ∂t L - S(t) dt ∂η [L - S(t)] Fave2Cpave2

[

{

Krw

∂2 (Pg1 - Pc) ∂ξ2

+

∂t

)

}

∂(Pg1 - Pc) ∂Krw ∂ξ ∂ξ

kave1 ξ dS ∂Pg1 Pg1 1 + × 2 T1 S (t) Fave1Cpave1 S(t) dt ∂ξ

{

∂2T1

}

2 M Cpg Κabsg ∂P g1∂ ln T1 Κ + rg ∂ξ ∂ξ ∂ξ2 2zR kave1 µg + FwCpw Κabs ∂(Pg1 - Pc) ∂T1 Κrw kave1 µw ∂ξ ∂ξ Κabs 1 × φ(1 - Sw) S2(t) Pg1 ∂2(Pg1 - Pc) Pg1∂Κrw ∂(Pg1 - Pc) K + µw rw µw ∂ξ ∂ξ ∂ξ2 2 2 2 ∂ P g1 ∂P 2g1 ∂ ln T1 1 1 1 ∂P g1 ∂Κrg + Κrg 2 + Κrg 2µg 2µg ∂ξ ∂ξ 2µg ∂ξ ∂ξ ∂ξ

{

+

(28)

-

}

}

FwCpw Κabsw ∂(Pg1 - Pc) ∂T1 Κ kave1 µw rw ∂ξ ∂ξ

(29)

(30)

Transformed equation for pressure in zone 2:

]

2 M Cpg Keff2 ∂Pg2 ∂ ln T2 + 2zR kave2 µg ∂η ∂η

∂η2 ∂2P2g2 ∂P2g2 ∂ 1n T2 Keff2 0.5 1 (31) ‚ ∂η ∂η φ(1 - Sh) µg [L - S(t)]2 ∂η2

[

]

Transformed equation for temperature in zone 2:

{

dS 1 ) dt ∆ΗFhφβ

ΚabsΚrw 1 ∂Pw PDM ΚabsΚrg 1 ∂Pg1 C - FwTDCpw µw S(t) ∂ξ zR pg µg S(t) ∂ξ PDM Κeff2 ∂Pg2 1 C + zR pg µg L - S(t) ∂η ∂T2 1 ∂T1 1 - kave1 + kave2 S(t) ∂ξ L - S(t) ∂η

}

(34)

Cpg ) heat capacity of gas (2227 J kg-1 K-1)27 Cpw ) heat capacity of water 4183 J kg-1 K-1)27 Cps ) heat capacity of sand (800 J kg-1 K-1)24 Cph ) heat capacity of hydrate (2200 J kg-1 K-1)15 kg ) heat conductivity of gas (0.03427 Wm-1 K-1)25 kw ) heat conductivity of water (0.6089 Wm-1 K-1)25 ks ) heat conductivity of sand (1.3 Wm-1 K-1)26 kh ) heat conductivity of hydrate (0.49 Wm-1 K-1)15 L ) length of simulation Pc ) capillary pressure Pg ) gas pressure Pw ) water pressure S ) location of the moving front Sg ) gas saturation Sw ) water saturation Sh ) hydrate saturation (0.19)4,7 T ) gas temperature (K) Ug ) superficial gas velocity Uw ) superficial water velocity Z ) compressibility of gas (0.88)21 ave ) average Greek Letters

∂Pg2 kave2 1 1 - η ∂Pg2 Pg2 ) + × 2 ∂t T2 [L - S (t)] Fave2Cpave2 L - S(t) ∂η

[

[

Symbol Definition

{

kave1 ∂T1 ∂2T1 ξ dS ∂T1 1 ) + + ∂t S(t) dt ∂ξ S2(t) Fave1Cpave1 ∂ξ2 2 M Cpg Κabsg ∂P g1 ∂ln T1 Κ + 2zR kave1 µg rg ∂ξ ∂ξ

+

dS ) dt

∂Pg2 ΚabsgΚrg 1 ∂Pg1 Κeff2 1 + µg L - S(t) ∂η µg S(t) ∂ξ (33) zRTD - (β - RD) φ βFh P DM

Nomenclature

Transformed equation for temperature in zone 1:

∂2T2

]

2 M Cpg Keff2∂Pg2∂ 1n T2 (32) 2zR kave2 µg ∂η ∂η

Energy balance at front:

Transformed equation for pressure in zone 1: ∂Pg1

∂η2

+

Mass balance at front:

The transformed governing equations using the coordinate transformation method are shown as follows. Here the subscript D means the variables at the dissociation moving front. Transformed equation for water saturation in zone 1:

∂Sw ξ dS ∂Sw 1 Kabs 1 × ) + ∂t φ µw S2(t) S(t) dt ∂ξ

∂2T2

]

φ ) porosity (0.2)4,7  ) mass fraction of methane per unit mass of hydrate (0.129)1 R ) water saturation in zone 1 (dissociated zone) β ) hydrate saturation in zone 2 (undissociated zone) Fg ) density of gas Fw ) density of water (1000 kg‚m-3) Fh ) density of hydrate (910 kg‚m-3)15 Fs ) density of sand (2200 kg‚m-3)24 Κabs ) absolute permeability of the porous media Κeffective ) effective permeability Κrg ) relative permeability of gas

2828

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Κrw ) relative permeability of water µg ) viscosity of methane gas (1.0 × 10-5 Pa‚s)27 µw ) viscosity of water (1.0 × 10 -3 Pa‚s)27 ∆H ) latent heat of hydrate dissociation23,24 Literature Cited (1) Sloan, E. D. Clathrate hydrates of natural gases; 2nd ed.; Marcel Dekker: New York, 1998. (2) Max, M. D.; Johnson, A. H.; Dillon, W. P. Economic Geology of Natural Gas Hydrate; Springer: Berlin, 2006. (3) Lee, S.-Y.; Holder, G. D. Methane Hydrates Potential as a Future Energy Source. Fuel Process. Technol. 2001, 71, 181-186. (4) 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. (5) 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. (6) Ahmadi, G.; Ji, C.; Smith, D. H. Natural gas production from hydrate dissociation: An axisymmetric model. J. Pet. Sci. Eng. 2007, 58, 245258. (7) Ji, C.; Ahmadi, G.; Smith, D. H. Natural gas production from hydrate decomposition by depressurization. Chem. Eng. Sci. 2001, 56, 5801-5814. (8) Ji, C.; Ahmadi, G.; Zhang, Wu.; Smith, D. H. Natural Gas Production From Hydrate Dissociation: A Comparison of Asymmetric Models. Proceedings of the Fourth International Conference on Gas Hydrates, Yokohama, Japan, 2002, May 19-23. (9) 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. (10) Moridis, G. J.; Collett, T. S.; Dallimore, S. R.; Satoh, T.; Hancock, S.; Weatherill, B. Numerical studies of gas production from several CH4 hydrate zones at the Mallik site, Mackenzie Delta, Canada. J. Pet. Sci. Eng. 2004, 43, 219-238. (11) Moridis, G. J.; Sloan, E. D. Gas production potential of disperse low-saturation hydrate accumulation in oceanic sediments. Energy ConVers. Manage. 2007, 48, 1834-1849. (12) Nazridoust, K.; Ahmadi, G. Computational modeling of methane hydrate dissociation in a sandstone core. Chem. Eng. Sci. 2007, 62, 61556177. (13) Sun, X.; Mohanty, K. K. Kinetic simulation of methane hydrate formation and dissociation in porous media. Chem. Eng. Sci. 2006, 61, 3476-3495.

(14) Selim, M. S.; Sloan, E. D. Heat and Mass Transfer During the Dissociation of Hydrate in Porous Media. AIChE J.1989, 35, 1049-1052. (15) Makogon, Y. F. Hydrates of Hydrocarbons; Penn Well Publishing Company: Tulsa, OK, 1997. (16) Bejan, A.; Kraus, A. D. Heat Transfer Handbook; John Wiley and Sons: New York, 2003. (17) Arastoopour, H.; Semrau, J. Mathematical Analysis of Two-Phase Flow in Low Permeability Porous Media. AIChE J. 1989, 35, 1710-1718. (18) Bear, J. Dynamics of Fluids in Porous Media; Dover Publications: Mineola, NY, 1988. (19) Corey, A. T. Mechanics of Immiscible Fluids in Porous Media; Water Resources Publication: Littleton, Colorado, 1994. (20) 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. (21) Smith, J. M.; Van Ness, H. C.; Abbott, M. M. Introduction to Chemical Engineering Thermodynamics, 5th ed.; McGraw-Hill: New York, 1996. (22) Masuda, Y.; Fujinaga, Y.; Naganawa, S.; Fujita, K.; Sato, K.; Hayashi, Y. Modelling and Experimental Studies on Dissociation of Methane Gas Hydrates in Berea Sandstone Cores. Proceedings of the 3rd International Conference on Gas Hydrates, Salt Lake City, UT, 1999. (23) Kang, S.-P.; Lee, H.; Ryu, B.-J. Enthalpies of dissociation of clathrate hydrates of carbon dioxide, nitrogen, (carbon dioxide + nitrogen), and (carbon dioxide + nitrogen + tetrahydrofuran). J. Chem. Thermodyn. 2001, 33, 513-521. (24) Lide, D. R. CRC Handbook of Chemistry and Physics, 85th ed.; CRC Press: Boca Raton, FL, 2005. (25) Wang, G. X.; Prasad, V.; Matthys, E. F. An Interface-Tracking Numerical Method for Rapid Planer Solidification of Binary Alloys with Application to Micro-segregation. Materials Sci. Eng. 1997, A225, 47-58. (26) Schiesser, W. E. The Numerical Method of Lines Integration of Partial Differential Equations; Academic Press: San Diego, CA, 1991. (27) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena, 2nd ed.; John Wiley & Sons: New York, 2002.

ReceiVed for reView October 16, 2007 ReVised manuscript receiVed January 30, 2008 Accepted February 12, 2008 IE071398B