Proceedings of the 7th International Conference on Gas Hydrates (ICGH 2011), Edinburgh, Scotland, United Kingdom, July 17-21, 2011.
MODELING DISSOCIATION AND REFORMATION OF METHANE AND CARBON DIOXIDE HYDRATE USING PHASE FIELD THEORY WITH IMPLICIT HYDRODYNAMICS Muhammad Qasim Department of Physics and Technology, Faculty of Mathematics and Natural Sciences University of Bergen Allégaten 55, 5007 Bergen NORWAY Bjørn Kvamme∗ Department of Physics and Technology, Faculty of Mathematics and Natural Sciences University of Bergen Allégaten 55, 5007 Bergen NORWAY Khuram Baig Department of Physics and Technology, Faculty of Mathematics and Natural Sciences University of Bergen Allégaten 55, 5007 Bergen NORWAY ABSTRACT Natural gas hydrates may dissociate or reform if they are exposed to phases which are under saturated with respect to hydrate with respect to water or hydrate formers. Contact with mineral surfaces will also lead to dissociation under the creation of a new interface bridging the hydrate and mineral. The kinetic rates of hydrate dissociation towards these undersaturated phase are important in the understanding of natural hydrates in sediments and the impact of contact with surrounding fluid phases and adsorbed phase (on mineral surfaces). Natural gas hydrates can also be produced though injection of CO2 because CO2 hydrate with CH4 left in small cavities will be more stable than the original natural gas hydrate. Common to this exchange process and the dissociations toward undereaturated phases is that the kinetic rates will depend on how fast the processes happen. If the dissociation rate per hydrate former is slower than self diffusion into the surrounding phase then the kinetic rate is not dependent on hydrodynamics. If - on the other hand - the dissociation is faster than the ability of surroundings to dilute the released molecules then bubbles or droplets will form and this will affect the phase transition rate and gravity forces needs to be included in the model. In this work an exchange of CO2 with CH4 hydrate is carried out using phase field theory considering variable density of the system coupled with hydrodynamics model. It is observed that the amount of the penetrated CO2 is about 60 % as that of the dissociated CH4. Keywords: Phase field theory, Natural gas hydrate, Hydrodynamics, Dissociation, Hydrate, exchange.
∗ Corresponding author: Phone: +4755583310 E-mail:
[email protected]
NOMENCLATURE Constant Constant Diffusion Coefficient [m2/s] Free energy functional [KJ] Free energy density [KJ/m3]
Gibbs free energy of inclusion of ∆ Component k in cavity type j [kJ]
Enthalpy [kJ] Cavity partition function of Component k in cavity type j & Concentration field mobility and phase field mobility [m3.s/J] Free energy scale [J.K/m3] Pressure [Pa] Pressure [Pa] Reference pressure [Pa] Molar gas constant [kJ/Kmol] Temperature [K]
Molar volume of the ith component [m3/mole] Average molar volume of pure water [m3/mole] Average molar volume of pure water [m3/mole] No. of type j cavities per water molecule Mole fraction Mole fraction of water Compressibility factor μ Chemical potential [KJ/mol] μ,# Chemical potential for water in empty hydrate structure [kJ/mol] $ Fractional occupancy of cavity k by component j % Activity Coefficient & Order parameter ' Molar density [Kg/m3] (
Number of type ) cavities per water molecule INTRODUCTION Gas hydrates are ice-like substances of water molecules encaging gas molecules (mostly methane) that form under high pressure and low temperature conditions within the upper hundred meters of the sub-seabed sediments [1]. These gas hydrates are widely distributed in sediments along continental margins, and harbor enormous
amounts of energy. Massive hydrates that outcrop the sea floor have been reported in the Gulf of Mexico [2]. Hydrate accumulations have also been found in the upper sediment layers of Hydrate ridge, off the coast of Oregon and a fishing trawler off Vancouver Island recently recovered a bulk of hydrate of approximately 1000kg [3]. Håkon Mosby Mud Volcano of Bear Island in the Barents Sea with hydrates is openly exposed at the ocean floor [4]. These are only few examples of the worldwide evidences of unstable hydrate occurrences that leaks methane to the oceans and eventually may be a source of methane increase in the atmosphere. Hydrates of methane are not thermodynamically stable at mineral surfaces. From a thermodynamic point of view the reason is simply that water structure on hydrate surfaces are not able to obtain optimal interactions with surfaces of calcite, quarts and other reservoir minerals. The impact of this is that hydrates are separated from the mineral surfaces by fluid channels. The sizes of these fluid channels are not known and are basically not even unique in the sense that it depends on the local fluxes of all fluids in addition to the surface thermodynamics. Stability of natural gas hydrate reservoirs therefore depends on sealing or trapping mechanisms similar to ordinary oil and gas reservoirs. Many hydrate reservoirs are in a dynamic state where hydrate is leaking from top by contact with groundwater/seawater which is under saturated with respect to methane. Dissociating hydrate degasses as bubbles if dissociation rate is faster than dilution in surrounding fluids and/or surrounding fluid is supersaturated. The kinetic rate depends on mass transport dynamics as well as thermodynamic driving force. Phase field theory is a power full tool to quantify this balance and provide a theoretical basis for development of simplified models for reservoir modeling tools. Gas hydrates are a potential source of energy as well as a potential solution for the reduction of CO2 emissions. CO2 hydrate is more stable than CH4 hydrate over substantial regions of pressure and temperature and a mixed hydrate of structure I in which CO2 fills large cavities and CH4 fills small cavities are stable over all regions of temperatures
and pressures. Gas hydrates have great capacity to store gases [5-7] and several investigations of potential for using hydrate phase for storage and transport have been conducted. The primary goal of this work is to develop a kinetic model for the replacement of CH4 with CO2 in natural gas hydrate reservoirs using modified Phase Field Theory [8]. This process is a favorable way to store a greenhouse gas (CO2) for long period of time and enables the ocean floor to remain stabilized even after recovering the methane gas [9]. THERMODYNAMICS Phase field theory model follows the formulation of Wheeler et al. [10], which historically has been mostly applied to descriptions of the isothermal phase transition between ideal binary-alloy liquid and solid phases of limited density differences. The hydrodynamics effects and variable density were incorporated in a three components phase field theory by Kvamme et al. [11] through implicit integration of Navier stokes equation following the approach of Qasim et al.[8]. The phase field parameter & is an order parameter describing the phase of the system as a function of spatial and time coordinates. The phase field parameter & is allowed to vary continuously from 0 to 1 on the range from solid to liquid. The solid state is represented by the hydrate and the liquid state represents fluid and aqueous phase. The solidification of hydrate is described in terms of the scalar phase field &* + , , , - where + , , , represents the molar fractions of CH4, CO2 and H2O respectively with obvious constraint on conservation of mass ∑ /+ 0 1. The field & is a structural order parameter assuming the values & 0 0 in the solid and & 0 1 in the liquid [12]. Intermediate values correspond to the interface between the two phases. The starting point of the three component phase field model is a free energy functional [11], : 89
F0
4 56 7 , *;&-, < ∑ , /+
: 8=>,?
@
'A ; B
; -, < CDE *&, + , , , , -F ,
(1)
which is an integration over the system volume, while the subscripts ), G represents the three components, ' is molar density depending on relative compositions, phase and flow. The bulk free energy density described as CDE 0 *&- < A1 B *&-IJ * + , , , , - < *&-K * + , , , , - .
(2)
The phase field parameter switches on and off the solid and liquid contributions J and K through the function *&- 0 & *10 B 15& < 6& , -, and note that *0- 0 0 and *1- 0 1. This function was derived from density functional theory studies of binary alloys and has been adopted also for our system of hydrate phase transitions. The binary alloys are normally treated as ideal solutions. The free energy densities of solid and liquid is given by OPQ
J 0 NJ '
(3)
,
K (4) K 0 NK ' . The thermodynamics for the hydrate system is discussed in more detail in the thermodynamics section below, with derivations of the free energies NJ for hydrate and NK for liquid state. Hydrate #PQ density ' is calculated using the formulation by Sloan et al. [13] K The liquid density ' for fluid phase is calculated as
KRED Q
0 S
/+
K,RED Q ' 0
+ UVWXY>Z T
(5)
,
where represents the molar volume of ith component. The molar volume is calculated using gas law 0
[ ] < \ ^ ][ _,`,
>a?
,
(6)
where represents the pressure and is compressibility factor calculated using SRK equation of state. For simplicity to avoid partial molar volume at infinite dilution the density of
liquid in aqueous phase is calculated as bcDdeDf
K
,
0 S <
/+
K,bcDdeDf '
0
1
(7)
bcDdeDf
K
,
where is the average molar volume of pure water. The function *&- 0 & , *1 B & , -⁄4 ensures a double well form of the CDE with a free energy scale 0 i1 B > l m < > n with jk
jk
*0- 0 *1- 0 0, where is the average molar volume of water. In order to derive a kinetic model we assume that the system evolves in time so that its total free energy decreases monotonically [12]. The usual equations of motion are supplemented with appropriate convection terms as explained by Tegze et al [14]. Given that the phase field is not a conserved quantity, the simplest form for the time evolution that ensures a minimization of the free energy is o op o> op
< *q. s-& 0 B *&, + , , , -
tR t
,
tR
< *q. s- 0 s. i *&, + , , , -s l t
,
where q is the velocity, 0 *1 B 0 i1 B
> l m jk
n jk
+ - u`
(8)
(9)
and
are the mobilities
associated with coarse-grained equation of motion which in turn are related to their microscopic counter parts. Where 0 J < *K B J -*&- is the diffusion coefficient. The details are given elsewhere [11]. An extended phase field model is formulated to account for the effect of fluid flow, density change and gravity. This is achieved by coupling the time evolution with the Navier Stokes Equations. The phase and concentration fields associates hydrodynamic equation as described by conti [1517]
'
wq oj op
ov op
0 B' ;. q,
< '*q. ;-q 0 'q < ;. .
(10)
Where q is the gravitational acceleration. ' is the #PQ density of the system in hydrate (' ) and liquid K ). Further (' 0 x < y.
(11)
is the generalization of stress tensor [15-18], x represents non-dissipative part and Π represents the dissipative part of the stress tensor. The expression for chemical potential of water in hydrate is {# 0 {,# B S ( |[ 71 < S F ,
(12)
This equation is derived from the macro canonical ensemble under the constraints of constant amount of water, corresponding to an empty lattice of the actual structure. Details of the derivation are given elsewhere [19] and will not be repeated here. {,# is the chemical potential for water in an empty hydrate structure and is the cavity partition function of component G in cavity type ). The first sum is over cavity types, and the second sum is over components G going into cavity type ). Here ( is the number of type ) cavities per water molecule. Fluid Thermodynamics The free energy of the fluid phase is assumed to have NK 0 ∑ /+ { RED Q ,
(13)
{RED Q 0 {∞ < |[* - ,
(14)
where { RED Q is the chemical potential of the ith component. The solubility of water is assumed to follow the Raoult’s law. The lower concentration of water in the fluid phase and its corresponding minor importance for the thermodynamics results in the following form of water chemical potential with some approximation of fugacity and activity coefficient:
where {∞ chemical potential of water at infinite dilution and is the mole fraction of water in the fluid phase. The chemical potential for the mixed fluid states is approximated as
{ RED Q
0
Ju},_D~d {
< |[* - ,
(15)
where ) represents CH4 or CO2. The details are available in Svandal et al. [20]. Aqueous Thermodynamics The free energy of the aqueous phase assumed as bcDdeDf
NK 0 ∑ /+ {
,
(16)
bcDdeDf
the chemical potential {
of aqueous phase has the general form derived from excess thermodynamics {# 0 { ∞ < |[* % ∞ - < * B -.
(17)
{ ∞ is the chemical potential of component ) in water at infinite dilution, % ∞ is the activity coefficient of component ) in the aqueous solution in the asymmetric convention (% ∞ approaches unity in the limit of becoming infinitely small). The chemical potentials at infinite dilution as a function of temperature are found by assuming equilibrium between fluid and aqueous bcDdeDf phasesA{RED Q 0 { I. This is done at low pressures where the solubility is very low, using experimental values for the solubility and extrapolating the chemical potential down to a corresponding value for zero concentration. The activity coefficient can be regressed by using the model for equilibrium to fit experimental solubility data. The chemical potential of water can be written as: _
{ 0 { < |[*1 B -% < * B -,
remain constant in the system. The values for temperature and pressure are taken at Nyegga cold seeps located on the edge of the Norwegian continental slope and the northern flank of the Storegga slide, on the border between two large oil/gas prone sedimentary basins – the Møre basin to the south and the Vøring basin to the north [2122]. The temperature and pressure condition is well inside the stability region of the guest molecules. The standard value of 9.8 m/s2 is assumed for gravity along the Y-axis of 2D geometry.
(18)
_
where { is pure water chemical potential. The strategy for calculating activity coefficient is given in [20]. RESULTS AND DISCUSSIONS The phase field model is implemented on a 2D geometry. This structure is used for the exchange of methane with carbon dioxide. Figure 1 show the circle of methane hydrate with blue color which is placed in the center surrounded by liquid carbon dioxide. The size of system is (500×500) grids with diameter of 200 grids for circular hydrate. One grid is equal to one angstrom (Å) and temperature (273.21 K) and pressure (63.90 bar)
Figure 1. Simulation at time zero, showing the initial picture of CH4 hydrate and liquid CO2 with 500x 500 grid points and a hydrate radius of 200 grid points, color codes represent φ = 0 and 1 for methane hydrate and CO2 liquid phases respectively. The simulation is run to 91.809E+06 total time steps this corresponds to the time of 91.809 ns Total number of grid points System area in m2 No. of time steps Total time in seconds CH4 mole fraction in hydrate CO2 mole fraction in Liquid
500×500 2.5E-15 91.809E+06 91.809E-09 0.14 1.0
Table 1. Simulation setup.
Figure 2. Structural phase parameter of the dissociating CH4 hydrate and exchange of CO2 at
different times. Between phases = 0 and 1 are interface values. Figure 2 shows that the methane initially starts to dissociate into the surrounding CO2 liquid. This is due to the driving force in terms of the change in chemical potential of methane in liquid phase and hydrate phase. CO2 is assumed to only enter the large cavities of structure I due to its size. The CO2 starts penetrating into the methane hydrate after some amount of methane has been released into the liquid phase. The hydrate size is reduced because of methane dissociation until 9.025 nano seconds and then it increased until 91.809 nano seconds due to reformation with CO2 penetration. This phenomenon is best explained by assistance of figure 3 which illustrates the reformation process of CO2 hydrate, where the kinetics of the liquid CO2 from its liquid phase transformation to solid phase at different stages is plotted. A thick interface between liquid and solid phase is highlighted with black circle in figure 3.
Figure 4. CO2 concentration inside hydrate as a function of radial distance at different times 0.0 ns, 6.4 ns, 10 ns, 22.5 ns, 40.0 ns and 77.84 ns. To observe the movement of methane from solid phase to liquid, the velocity on the interface is determined by tracking the & values which is used to calculate the dissociation rate until 91.809 ns using the following equation [3]: 0
vZ
i
l.
.+#
(19)
where is the dissociation rate (moles/m s), is hydrate radius shrinkage rate (m/s), '#PQ is 2
density of hydrate (kg/m ), is molar weight of the guest (kg/moles) and is Hydrate number. The calculated flux profile is plotted in figure 5. The initial value of flux was high due to the initial relaxation of a system into a physically realistic interface. The actual dependency of formation rate on driving forces is illustrated in figure 5. 3
Figure 3. CO2 concentration as a function of distance and times 0 ns, 6.4 ns, 10 ns, 22.5 ns, 40.0 ns and 77.84 ns. The concentration of CO2 in the hydrate as a function of radial distance is shown in figure 4. The totally converted liquid CO2 is shown by point A on right side of figure 4. This point was used to make an estimate of converted CO2 which was approximated about 0.12 mole fraction. The system formed an interface with liquid phase between A to B resulting into an initial relaxation of the system into a physically realistic interface.
Figure 5. Carbon dioxide flux as a function of time. The rate is decreasing gradually after this relaxed point on the curve. One reason for this is the decrease in thermodynamic driving force, which is proportional to the increasing chemical potential inside the hydrate, which is still filled with methane. The noise seen on the calculated curves is due to grid resolution effects. The interface in this simulation perfectly follows a power law
which is proportional to square root of time (α t1/2), which according to Fick's law indicates a diffusion controlled process. The total number of CO2 molecules penetrated into methane hydrate is 1.5308E+03. The simulated results have been extrapolated to 3.1536E+20 nano seconds from 91.809 nano seconds according to this observation and the maximum time in figure 6 is equal to 104 years. This is done to see the possible trend on a long time scale.
liquid CO2 is increasing. Higher concentration of methane was observed to have diffused into the liquid phase by forming a transition zone at the end of the simulation. Some methane molecules have diffused rapidly into the liquid phase because some of the vertices of water cages in the hydrate has dissociated shown on the right side of figure 8.
Figure 8. Methane concentration inside the hydrate as a function of radial distance at different times. Figure 6. Extrapolation of carbon dioxide flux up to 10000 years. Due to the length of the time scale the values were plotted in the figure with 100 years of time intervals. After 10000 years the reformation rate of CO2 is 2.198E-12 moles/m2s which corresponds to 6.932E-05 moles/m2yr.
The interface became stable after some time, as illustrated in fig. 8. This can be explained through the definition of solubility which is a measure of how much solute (methane) will dissolve in all the solvent (carbon dioxide). The methane will not completely dissolve in carbon dioxide because of the polar solvents molecules separate the molecules of other polar substance. There is no thermodynamic equilibrium between the wall of the hydrate and that of the fluid in sharp stable interface. This sharp interface is shown in figure 9 which illustrate the concentration of methane in liquid and solid phases.
Figure 7. Methane concentration as function of length. Figure 7 illustrate the change of methane concentration with respect to position in simulation box while methane diffuses into liquid phase and CO2 penetrating into the large cavities of methane hydrate. At time t=0 nano seconds there is no change in methane concentration but as the time passes to t=0.4 interface between liquid and solid phase is developed which represents the presence of methane in liquid and solid phases. The interface thickness is increasing as the time is increasing, which is highlighted with black circles in figure 7, and also the methane concentration in
Figure 9. Methane concentration as a function of length at different times. Figure 10 shows the diffusion of methane from hydrate phase to liquid phase. This is sampled by tracing the movement of methane concentration from the interface between liquid and solid.
Figure 10. Methane dissociation rate. Initially the value of methane flux is higher because of initial relaxation of the system as discussed in the case of carbon dioxide flux. To show the actual dependency of dissociation on driving forces the close look on the curve is shown in figure 10. The flux is gradually decreasing after this relaxation is because of the decrease in thermodynamic driving forces which is proportional to the increasing methane concentration in the surrounding liquid. The value of methane at the end of simulation is 9.297 moles/m2s. This high value is because system is not reached at equilibrium. The total number of 2.5545E+03 methane molecules dissociated after 91.809. The interface allows almost perfectly the power law which is proportional to square root of time (α t1/2) showing a diffusion control process.
Figure 12. Methane and carbon dioxide concentration profile in solid and liquid phase. The CH4 and CO2 profiles are going opposite directions which is due to dissociation of methane from hydrate to liquid phase and CO2 reformation which refilling the large cavities of the methane hydrate. This is best illustrating by figure 13. Over the period of simulations, the comparison of CH4 dissociation and CO2 reformation resulted in approximately 60 percent of methane exchanged while CO2 is stored as CO2 hydrate. Theoretically, the maximum exchange could be 75 % corresponding to CO2 filling all the large cavities.
The simulation is extrapolated to 3.1536E+20 nano seconds from 91.809 nano seconds which is equal to 10000 years. The flux of the methane after 10000 years is 2.954E-14 moles/m2s which corresponds to 9.316E-07 moles/m2yr.
Figure 13. Methane and carbon dioxide concentration profile inside the hydrate at different times.
Figure 11. Extrapolation of methane flux up to 10000 years. The comparison between CH4 and CO2 is shown in figure 12 where CH4 and CO2 concentration profile were combined to illustrate different phases achieved during the exchange process.
CONCLUSION Phase field simulation with more appropriate description of density dependencies in the phase field part as well as in the hydrodynamic parts [8] has been applied to model the exchange of CH4 with CO2 from natural gas hydrate at conditions corresponding to hydrates in Nyegga. The kinetic data attained are examples of important results which will be useful in the modeling and optimization for the production of methane from
hydrate reservoir as well as sequestration of CO2. As expected it was observed that the mole fraction of CO2 in the hydrate phase increases, while that of CH4 decreases with increasing time. Within the limited simulations time approximately 60% of initial CH4 in hydrate was exchanged. REFERENCES [1] Kvenvolden K. A, Rogers B. W. Gaia’s breath—global methane exhalations. Marine and Petroleum Geology 2005;Vol.22: 579-590. [2] MacDonald I. R, Guinasso N. L, Sassen Jr. R, Brooks J. M, Lee L, Scott K. T. Gas hydrate that breaches the sea floor on the continental slope of the Gulf of Mexico. 1994:699-702. [3] Rehder G, Kirby S. H, Durham W. B, Stern L. A, Peltzer E. T, Pinkston J, Brewer P. G. Dissolution rates of pure methane hydrate and carbon-dioxide hydrate in undersaturated seawater at 1000-m depth. Geochimica et Cosmochimica Acta 2004; 68(2):285-292. [4] Egorov A. V, Crane K, Vogt P. R, Rozhkov A. N, Shirshov P. P. Gas hydrates that outcrop on the sea floor: stability models. Geo-Marine Letters 1999;Vol.19:68-75. [5] Saji A, Yoshida H, et al. Fixation of carbondioxide by clathrate-hydrate. Energy Conversion and Management 1992;33(5-8):643-649. [6] Yamasaki A, Teng H, et al. CO2 hydrate formation in various hydrodynamics conditions. Gas Hydrates: Challenges for the Future 2002; 912:235-245. [7] Lee J. W, Chun M. K, et al. Phase equilibria and kinetic behavior of CO2 hydrate in electrolyte and porous media solutions: Application to ocean sequestration CO2. Korean Journal of Chemical Engineering 2002; 19(4):673-678. [8] Qasim M, Kvamme B, Baig K. Phase field theory modeling of CH4/CO2 gas hydrates in gravity fields. International Journal of Geology 2011; 5(2):48-52. [9] Caldeira K, Wickett M. E. Ocean model predictions of chemistry changes from carbon dioxide emissions to the atmosphere and ocean. J. Geophys. Res. 2005; 110:12. [10] Wheeler A. A, Boethinger W. J, McFadden G. B, Phase field model for isothermal phase transitions in binary alloys. Physical Review A 1992; Vol. 45:7424-7439. [11] Kvamme B, Svandal A, Buanes T, KuznetsovaT. Phase field approaches to the kinetic modeling of hydrate phase transitions. AAPG Memoir 2009 ;Volume 89:758-769
[12] Svandal A. Modeling hydrate phase transitions using mean field approaches. University of Bergen 2006:1-37. [13] Sloan E. D, Koh C. A. Clathrate hydrates of natural gases (3rd ed.). Boca Raton, FL: CRC Press, 2008. [14] Tegze G, Gránásy L. Phase field simulation of liquid phase separation with fluid flow. Material science and engineering 2005; vol.413-414:418422. [15] Conti M. Density change effects on crystal growth from the melt. Physical Review 2001; vol.64:051601. [16] Conti M, Fermani M. Interface dynamics and solute trapping in alloy solidification with density change. Physical Review 2003; vol.67:026117. [17] Conti M. Advection flow effects in the growth of a free dendrite. Physical Review 2004; vol.69:022601. [18] Baig K, Qasim M, Kivelä P. H, Kvamme B. Phase field theory modeling of methane fluxes from exposed natural gas hydrate reservoirs. American Institute of Physics 2010, in press. [19] Kvamme B, Tanaka H. Thermodynamic stability of hydrates for ethane, ethylene and carbon dioxide. J. Chem. Phys. 1995; vol. 99:7114-7119. [20] Svandal A, Kuznetsova T, Kvamme B. Thermodynamic properties and phase transitions in the H2O/CO2/CH4 system. Fluid Phase Equilibria 2006; vol.246:177-184. [21] Bunz S, Mienert J, Berndt C. Geological controls on the storegga gas-hydrate system of the mid-Norwegian continental margin. Earth and Planetary Science Letters 2003; 209(3-4):291-307. [22] Chen Y, Haflidason H, Knies J. Methane fluxes from pockmark area in Nyegga, Norwegian Sea. In: International Geological Conference, August 2008.
