Kinetics of solid hydrate formation by carbon dioxide: Phase field theory of hydrate nucleation and magnetic resonance imagingy
B. Kvamme,aA. Graue,aE. Aspenes,aT. Kuznetsova,aL. Gra´na´sy,bG. To´th,bT. Pusztaiband G. Tegzeb
a Department of Physics, University of Bergen, Alle´gt. 55, 5007 Bergen, Norway
b Research Institute for Solid State Physics and Optics, H-1525 Budapest, POB 49, Hungary
Received 15th September 2003, Accepted 14th November 2003 F|rst published as an Advance Article on the web 8th December 2003
In the course of developing a general kinetic model of hydrate formation/reaction that can be used to establish/
optimize technologies for the exploitation of hydrate reservoirs, two aspects of CO2hydrate formation have been studied. (i) We developed a phase field theory for describing the nucleation of CO2hydrate in aqueous solutions. The accuracy of the model has been demonstrated on the hard-sphere model system, for which all information needed to calculate the height of the nucleation barrier is known accurately. It has been shown that the phase field theory is considerably more accurate than the sharp-interface droplet model of the classical nucleation theory. Starting from realistic estimates for the thermodynamic and interfacial properties, we have shown that under typical conditions of CO2formation, the size of the critical fluctuations (nuclei) is comparable to the interface thickness, implying that the droplet model should be rather inaccurate. Indeed the phase field theory predicts considerably smaller height for the nucleation barrier than the classical approach. (ii) In order to provide accurate transformation rates to test the kinetic model under development, we applied magnetic resonance imaging to monitor hydrate phase transitions in porous media under realistic conditions. The mechanism of natural gas hydrate conversion to CO2-hydrate implies storage potential for CO2in natural gas hydrate reservoirs, with the additional benefit of methane production. We present the transformation rates for the relevant processes (hydrate formation, dissociation and recovery).
I. Introduction
The worldwide amounts of carbon bound in gas hydrates are conservatively estimated to total twice the amount of carbon to be found in all known fossil fuels on Earth.1 Assuming the lattice is filled to the maximum capacity, dissociation of 1 m3 gas hydrate may release about 164 m3 methane under standard temperature and pressure (STP) condition.2
CO2hydrate is significantly more stable thermodynamically than methane hydrate. Storage of CO2in hydrate reservoirs through replacement of natural gas is therefore considered as an interesting option for safe long terms storage of CO2, which can be economically feasible due to the produced natural gas.
The efficiency of any exploitation strategy based on CO2 depends on the composite dynamics of the system, where knowledge of the kinetics of hydrate reformation is crucial.
Storage of CO2in aquifers is another option for reducing CO2 emissions to the atmosphere. This option is already in use outside the coast of Norway, where CO2from the Sleipner field is being injected into the Utsira formation. Similar storage option may be used outside the northern part of Norway. In these regions the seabed temperature may reach as low as 1C and there may be regions of the reservoirs beneath that are inside the hydrate stability zones. Therefore, hydrate may form homogeneously from dissolved CO2or at the interface between free CO2and water.
Development/optimization of the related technologies requires a detailed understanding of the processes involved, a knowledge that can be gained by combined experimenting and modeling. Two essential ingredients of such an approach are the collection of reliable experimental information on the transition rates, and the development of appropriate
theoretical tools that can describe all stages of the processes involved. For example, one of the least understood steps of hydrate formation is the early nucleation stage, in which nan- ometer sized hydrate crystallites formviathermal fluctuations.
Accordingly, in this paper, we apply two modern approaches to study the hydrate phase transitions involved in the storage of CO2in reservoirs. First, we present a model for the nucleation of CO2 hydrate under conditions specific to such reservoirs. In this we rely on the phase field theory (PFT); one of the most potent methods developed recently to model solidification in binary, ternary and multi-component melts. Over the past decade it has demonstrated its ability to describe complex solidification morphologies3 including ther- mal and solutal dendrites4–7and eutectic/peritectic fronts.8–12 In addition, we demonstrate the accuracy of our approach on the hard-sphere system, for which all essential properties are well known, making it an ideal testing ground for theory.
Secondly, we apply magnetic resonance imaging13(MRI) to monitor hydrate phase transitions. We present experimental results for the formation and dissociation of CO2hydrate in real porous media.
II. Nucleation theory
The freezing of homogeneous undercooled liquids starts with the formation of heterophase fluctuations whose central part shows crystal-like atomic arrangement. Those heterophase fluctuations that exceed a critical size (determined by the inter- play of the interfacial and volumetric contributions to the clus- ter free energy) have a good chance to reach macroscopic dimensions, while the smaller ones decay with a high probabil- ity. (Heterophase fluctuations of the critical size are called nuclei.) The description of the near-critical fluctuations is pro- blematic even in single component systems. The main difficulty y Presented at the 3rd International Workshop on Global Phase
Diagrams, Odessa, Ukraine, September 14–19, 2003.
PCCP
www.rsc.org/pccp
R E S E A R C H P A P E R
DOI:10.1039/b311202k
is that the typical size of the critical fluctuations, forming on the human time scale, is about 1 nm, comparable to the thick- ness of the crystal-liquid interface, which in turn extends to several molecular layers.14 Therefore, the droplet model of the classical nucleation theory, which relies on a sharp inter- face and bulk crystal properties, is inappropriate for such fluc- tuations. Field theoretic models, that predict a diffuse interface, offer a natural way to handle such a situation.15 For example, in recent works, the phase field theory has been shown to describe such fluctuations quantitatively.16,17
II.1 Phase field theory of nuclei
Our starting point is the standard phase field theory of binary alloys as developed by several authors.3,18 In the present approach, the local state of the matter is characterized by two fields; a structural order parameter, m, called the phase field, that describes the transition between the disordered liquid and ordered crystalline structures, and a conserved field, w, which may be either the coarse-grained density, r, or the solute concentration,c.
The structural order parameter can be viewed as the Fourier amplitude of the dominant density wave of the time averaged singlet density in the solid. As pointed out by Shen and Oxtoby19if the density peaks in the solid can be well approxi- mated by Gaussians placed at the atomic sites, all Fourier amplitudes can be expressed uniquely in terms of the ampli- tude of the dominant wave, thus a single structural order para- meter suffices. Here we takem¼0 in the solid andm¼1 in the liquid. We assume mass conservation, which implies that the integral of the conservative fields over volume is a constant.
The Helmholtz free energy of the system is a functional of these fields:
F ¼ Z
d3rf12e2TðHmÞ2þfðm;wÞg ð1Þ
where e is a constant, T is the temperature, and f(m,w) is the local free energy density. The first term on the right hand side is responsible for the appearance of the diffuse interface. The local free energy density has the form f(m,w)¼wTg(m)þ[1p(m)] fS(w)þp(m) fL(w), where the
‘‘ double well ’’ and ‘‘ interpolation ’’ functions have the forms g(m)¼14m2(1m)2and p(m)¼m3(1015mþ6m2), respec- tively, that emerge from the thermodynamically consistent for- mulation of the PFT,20 w is the free energy scale, while the Helmholtz free energy densities of the homogeneous solid and liquid,fSandfL, depend on the local value of the conser- vative field. These relationships result in a free energy surface that has two minima, whose relative depth depends on the deviation from equilibrium.
Being in unstable equilibrium, the critical fluctuation (the nucleus) can be found as an extremum of this free energy func- tional,18–20 subject to the solute conservation constraint dis- cussed above. To impose this constraint one adds the volume integral over the conserved field times a Lagrange multiplier, l, to the free energy:lR
d3rw(r). The field distributions, that extremize the free energy, have to obey the appropriate Euler–Lagrange (EL) equations, which in the case of such a local functional take the form
dF dw¼dc
dw¼ Hdc
dHw¼0; ð2Þ
where dF/dwstands for the first functional derivative of the Helmholtz free energy with respect to the field w, while cis the total free energy density. The EL equations have to be solved assuming that unperturbed liquid exists in the far field, while, for symmetry reasons, zero field gradients exist at the
center of the fluctuations. Under such conditions, the Lagrange multiplier can be identified asl¼ (@c/@w)r!1.
Assuming spherical symmetry, which is reasonable consider- ing the low anisotropy of the crystal– liquid interface at small undercoolings, the EL equations take the following form:
e2T d2m dr2 þ2
r dm
dr
( )
¼wTg0ðmÞ þp0ðmÞffLfSg; ð3aÞ
and
0¼wTgðmÞ þ½1pðmÞdfS
dwþpðmÞdfL dwdfL
dw
r!1
:
ð3bÞ Here0stands for differentiation with respect to the argument of the function. The last term in eqn. (3b) originates from the Lagrange multiplier. Since the right hand side of eqn. (3b) is a function of fieldswandm, it provides the implicit relation- shipw¼w(m). Accordingly, eqn. (3a) is an ordinary differen- tial equation for m(r). This equation has been solved here numerically using a fourth order Runge–Kutta method. Since mand dm/drare fixed at different locations, the central value of m that satisfies m!m1¼1 for r! 1, has been deter- mined iteratively. Having determined the solutions m(r) and w(r), the work of formation of the nucleusW* can be obtained by inserting the solution into the free energy functional. Pro- vided that the bulk free energy densities,fS(w) and fL(w), are known, the only model parameters, we need to fix to evaluate W*, arewande. These model parameters are related to the interface thickness and the interfacial free energy,6,21thus these quantities can be used to determinewande, and calculateW*
without adjustable parameters.
The steady state nucleation rate (number of nuclei formed in unit volume and time),JSS, can be calculated as
JSS¼J0expfW=kTg; ð4Þ using the classical nucleation prefactor,22 J0, verified experi- mentally on oxide glasses.23
II.2 Classical droplet model (CDM)
For comparison with our non-classical approach, we calculate the height of the nucleation barrier using the sharp interface droplet model of the classical nucleation theory. In this approach, the free energy of heterophase fluctuations of radius R is given as WCDM¼ (4p/3)R3Dgþ4pR2G1, where Dg¼r(mLmS) is the volumetric Gibbs free energy difference between the bulk liquid and solid,mLandmSare the chemical potential of the liquid and the solid, and G1 is the free energy of the equilibrium planar interface. Then, the free energy of the critical fluctuation of radiusR*CDM¼2G1/Dg isWCDM¼(16p/3)G31/Dg2. This approach is expected to be accurate if only the interface thickness is negligible with respect to the size of the critical fluctuation.
II.3 Material properties
We apply these models for crystal nucleation in the hard- sphere system, and for CO2 hydrate nucleation in aqueous CO2solutions.
In the case of the hard-sphere system, the virtually exact equations of state for the solid and liquid phases24 deduced from molecular dynamics (MD) simulations have been inte- grated to obtainfS(r) andfL(r). Owing to the complex form of these functions, eqn. (3b) could only be inverted numeri- cally. The interface thickness for the phase field was evaluated from the cross-interfacial variation of the height of the singlet density peaks determined by MD simulations.25,26The model parameters were fixed in equilibrium so that the free energy27
and thickness of the (111), (110), and (100) interfaces from molecular dynamics were recovered. Considering that due to the small anisotropy the equilibrium shape is expected to be nearly spherical, we also performed calculations with an aver- age of the interface free energy over the known orientations.
In the case of hydrate nucleation, the molar Gibbs free energy of the aqueous CO2 solution has been calculated as GL¼(1c)GL,WþcGL,CO2, where cis the mole fraction of CO2. The partial molar Gibbs free energy of water in solution has been obtained as GL,W¼G0L;WþRTln[(1c)gL,W(c)], where the free energy of pure water has been calculated as
G0L;W¼X3
i¼0
ki
Ti; ð5Þ
with coefficientskitaken from ref. 28. HereRis the gas con- stant and gL,W is the activity coefficient of water in solution.
The partial molar free energy of CO2 in solution is GL,CO2¼G1L;CO2þRTln[cgL,CO2(c)], where the molar free energy of CO2 at infinite dilution, G1L;CO2¼ 19.67 kJ/mol, has been taken from molecular dynamics simulations.29 The temperature and pressure dependent activity coefficient of CO2in aqueous solution deduced from CO2solubility experi- ments of Stewart and Munjal30have been fitted using the form
lngL;CO2¼X5
i¼0aiðTÞ½lnxi; ð6Þ withai(T) given by third order polynomials. The activity coef- ficient of water, gL,W in aqueous solution has been obtained from eqn. (6)viathe Gibbs–Duhem relationship.
The Gibbs free energy of the hydrate is given by GS¼(1c)GS,WþcGS,CO2. Owing to the lack of experimen- tal information, the partial molar quantities have been calcu- lated using the model described in ref. 28. For water and CO2we use the relationshipsGS,W¼G0S;WþRT(3/23)ln(1y), ), andGS,CO2¼GincS;COþRTln[y/(1y)], respectively, where the hole occupancy isy¼c/(3/23). Here, the partial molar Gibbs free energies of the empty clathrate, G0S;W, and that of guest inclusion, GincS;CO, are given by eqn. (6), with the appropriate kitaken from ref. 28.
Following other authors,31we approximate the free energy of the hydrate–solution interface by that of the ice–water inter- face, taken from the work of Hardy, 29.10.8 mJ m2.32 Owing to a lack of information on the CO2hydrate/aqueous solution interface, we use the 10%–90% interface thickness,d, (the distance on which the phase field changes between 0.1 and 0.9) as an adjustable parameter in the calculations. Mole- cular dynamics simulations on other clathrate hydrates indicate that the full interface thickness is about 2–3 nm,33that corres- ponds to roughlyd1–1.5 nm. Assuming that similar values apply for the CO2hydrate, we varydin the 0.5–1.5 nm range.
The computations were performed under conditions typical for the seabed reservoirs,i.e. T¼274 K,p¼15 MPa (1500 m depth), furthermore, we assume that water has been satu- rated by CO2(c¼0.033, obtained by extrapolating the rele- vant data by Teng and Yamasaki34). These experimental data are for synthetic average seawater. The salinity of ground- water in reservoirs may vary from close to zero up to seawater salinity in regions where the penetration of seawater dominates the salinity.
III. Experimental
Monitoring of hydrate phase transitions in real porous media (a sandstone core plug) has been performed using Magnetic Resonance Imaging.13The experimental setup is shown sche- matically in Fig. 1. The core length is 9.4 cm and the core dia- meter is 1.5 cm. Porosity was measured at 22.1% and the brine permeability at 1050 mD.
The core was saturated with brine, placed inside the core holder in the MRI tomograph and pressurized applying 80 bars confinement pressure. In an attempt to reproduce condi- tions that are, or would be, present in nature, i.e. gas with water present, a gas-water displacement process was initiated.
A liquid gas phase, CO2at 63 bars and 25C, was introduced at one end of the 100% brine saturated core sample. The liquid gas-water displacement was imaged. The core holder was then cooled to 2.0C and hydrate began to form. Imaging the change in the free water saturation as function time provides an indirect measure of the hydrate formation and disassocia- tion during temperature variations.
IV. Results and discussion
IV.1 Nucleation
A. Hard-sphere system.Fixing the model parameterseand wso that the free energy and the 10%–90% thickness of the equilibrium interfaces are recovered, the grand potential sur- face (Fig. 2), and the phase field and coarse-grained density profiles for the equilibrium solid-liquid interface have been determined. The profiles predicted for the (111) interface are Fig. 1 The core plug inside the MRI tomograph (left) has been cooled using a coolant fluid (right) that controls the temperature with an accuracy of 0.1C. The arrows show the flow of the fluids (coolant, confinement fluid for the pressure vessel and CO2).
Fig. 2 Grand potential density surface (o¼fmLr) for the hard- sphere system in equilibrium, evaluated with the free energy scalew corresponding to the (111) interface. (m–phase field,f–volume frac- tion.) Note, that the minima at (f¼0.494,m¼1) and (f¼0.546, m¼0) correspond to the equilibrium liquid and solid (fcc) phases.
The calculation has been performed for a colloidal suspension of hard-sphere diameters¼890 nm, a particle size comparable to that in many experiments. The small value of the grand potential density follows from the large ‘‘ molecular ’’ diameter. (T¼293 K.)
compared with their counterparts from MD simulations in Figs. 3. For the sake of this comparison, all quantities are pre- sented in a normalized form, (XXmin)/(XmaxXmin), where Xis either the density peak height, filtered density, 1m, or the coarse-grained density, whileXminandXmaxare the mini- mum and maximum values of these quantities. For the (111) and (100) interfaces, and to a smaller extent for the (110) inter- face, the structural order parameter profiles are in a good agreement with the cross-interfacial variation of the density peak height from MD. Note that while dhas been fitted, the agreement depends also on the shape/symmetry of the pre- dicted profile. It appears that the nearly symmetric structural order parameter profile of the PFT approximates reasonably well the behavior seen in the simulations. Remarkably, the density functional theories by Shen and Oxtoby,19 and Gra´na´sy and Pusztai35derived for the fcc structure predict sig- nificantly more asymmetric structural order parameter profiles for the crystal-liquid interfaces, which are sharper on the crys- tal side and have a long tail on the liquid side. Therefore, their match to the simulated profiles is apparently less satisfactory.
It is an intriguing question why the quartic form ofg(m), that is consistent with the symmetries of the base centered cubic (bcc) structure,36leads to a better fit to the simulation profile then those obtained by using free energy functionals consistent with the fcc symmetries. A distinct possibility is the presence of a bcc-like layer at the interface as reported in the case of the Lennard-Jones system.37
Summarizing, with model parameters fixed as described above, the PFT reproduces exactly the thermodynamics of the bulk phases, the interface free energy, and quite reasonably both the structural order parameter profile and the density profile. Therefore, we may calculate the height of the nuclea- tion barrier with some confidence.
The reduced nucleation barrier heights,W*/kT, calculated with eandw that fit to the properties of the (111) and (110) interfaces and to the average interface properties are shown in Fig. 4 for the phase field theory. For comparison, the predic- tions of the classical nucleation theory and the exact results from Monte Carlo simulations38 are also presented. The W*/kT vs. initial liquid volume fraction curves predicted by the PFT using the properties of the (111) and (110) interfaces envelope the MC simulations, while the predictions with the average interface properties fall close to the MC results. (The initial volume fraction of the liquid is defined as f1¼ps3r1/6, wherer1is the initial density of the liquid.)
The radial order parameter and density profiles for the criti- cal fluctuations (nuclei) corresponding to the volume fractions used in MC simulations of ref. 38 are shown in Fig. 5. Where available (f1¼0.5207) these profiles are in a reasonable agreement with the size of the critical fluctuation from MC simulation (cf. our Fig. 5 and Fig. 3 of ref. 38). Note the lack of bulk crystalline properties even at the centers of the critical fluctuations. As a result, the classical droplet model signifi- cantly underestimates the height of the nucleation barrier for all orientations. Since the nucleation rate is proportional to exp(W*/kT), these differences amount to orders of magni- tude. For example, the CDM calculation made with the aver- age interface free energy overestimates the nucleation rate by three to five orders of magnitude.
These findings for the hard-sphere model system indicate that our non-classical approach is considerably more accurate
Fig. 3 Reduced interfacial profiles for the (111) interface. (Solid lines: predictions of the phase field theory; heavy solid line: structural order parameter, 1m; light solid line: density. Dashed lines: results from molecular dynamics;25,26 heavy dashed line: height of singlet density peaks; light dashed line: filtered density profile. Note that the heavy lines should be compared to each other, as the light ones.)
Fig. 4 Reduced nucleation barrier heightvs.volume fraction of the initial liquid. [Upper three lines (solid): phase field theory (PFT); lower three lines (dashed/dotted): classical droplet model–CDM.] Within these triplets of lines for PFT and CDM, the upper and lower curves were calculated using the physical properties of the (110) and (111) interfaces, respectively (as indicated by caption in the figure). The pre- diction for the (100) interface (not shown here) falls slightly below the results obtained with the average interface properties (heavy lines). The results of Monte Carlo simulations38are denoted by triangles.
Fig. 5 Radial 1mprofiles (solid lines) and reduced volume fraction profiles (dashed lines) for critical fluctuations that correspond to the initial liquid volume fractions f1¼0.5207, 0.5277, and 0.5343.
The PFT calculations were performed witheandwcorresponding to the average interface properties (G1s2/kT¼0.6133 andd/s¼5.410).
HerefS(¼0.5796, 0.5875, and 0.5948, respectively) is the volume fraction of the solid that is in mechanical equilibrium with the initial liquid (the crystal that has the same pressure as the liquid of volume fraction f1).
than the classical droplet model, a conclusion supported by earlier evidence on the Lennard-Jones and ice–water systems.16,17
B. Formation of CO2hydrate in aqueous solution.The free energy surface, corresponding to e andw that reproduce the ice-water interface free energy and a 10%–90% thickness of 1 nm at T¼274 K, is shown in Fig. 6. The radial structural order parameter and concentration profiles of the critical fluc- tuations forming in saturated aqueous CO2 solution (c¼0.033) at the same temperature are shown in Fig. 7 as a function of the interface thickness,d. Remarkably, at the rea- listic interface thickness (d¼1.0 to 1.5 nm) the bulk crystalline structure is not yet established even at the center of the nucleus, indicating that the nucleus is softer (the molecules have larger amplitude of oscillations around the crystal sites) than the bulk crystal. Despite these, we have almost full hole-occupancy in the central part, c¼0.1235 or y¼0.946.
Furthermore, the interface thickness for the concentration field is far sharper than for structure. It extends to only a few A˚ , which is consistent with the picture that the nucleus is a small piece of hydrate crystal embedded into the solution, however, built of somewhat distorted H2O cages seen at the interface in MD simulations with realistic potentials39 (Fig. 8). It is remarkable, that the interface for the solute falls close to the classical radiusRCDM¼1.76 nm.
In agreement with previous studies on other systems,16,17a substantial difference can be seen between predictions made for the height of the nucleation barrier by the phase field the- ory and the classical droplet model (Table 1). This difference emerges from the fact that the interface thickness is compar- able to the size of the critical fluctuation. As one should expect, for decreasing interface thickness, the PFT prediction for the height of the nucleation barrier converges to that by the sharp interface CDM (WPFT!WCDMford!0; see Fig. 9), indicat- ing the coherency of our results. These findings imply that, similarly to many other systems, in hydrate forming systems the classical droplet model of crystal nuclei is rather inaccurate and should be replaced by more advanced approaches such as the phase field theory.
We note finally, that the predictedW* is rather sensitive to the values ofG1anddused to fix the model parameterseand w. Thus, accurate results may only be expected if these input properties are known with a high accuracy. Unfortunately, none of these input properties is available for the CO2
hydrate/solution interface (the values used here and in other
work31are only rough estimates). Therefore, further effort is needed to establish accurate nucleation rates.
The experimental determination of the input properties is susceptible to errors. In the case of transparent substances, meticulous experiments using the grain boundary groove tech- nique may reduce the error of the interfacial free energy to 3%.32In ideal case, a comparable error can be established for the thermodynamic properties. These add up to 15%
error inW*, corresponding to about 2.5 orders of magnitude in the nucleation rate atW*/kT¼40. Another possibility is to evaluate the interfacial and thermodynamic properties from MD simulations25,27,40 with realistic interaction potentials.
Then, in principle, virtually exact thermodynamic properties can be used (for example, for the hard sphere fluid the first eight virial coefficients have been determined), and the inter- face free energy can be evaluated with an error as low as 1.52%.27 Accordingly, the overall uncertainty of W* is 56%, which corresponds to about an order of magnitude uncertainty in the nucleation rate at W*/kT¼40. This appears to be the maximum precision one may achieve for Fig. 6 Helmoltz free energy density surface for the CO2hydrate sys-
tem at 274 K, evaluated with free energy scalewcorrespondingd¼1 nm andG1¼29.10.8 mJ m2. (m: phase field,c: mole fraction of CO2.)
Fig. 7 Radial 1mprofiles (solid lines) and CO2concentration pro- files (dashed lines) for critical fluctuations calculated atT¼274 K, for 10%–90% interface thicknessd¼(a) 0.5, (b) 1.0 and (c) 1.5 nm. For comparison the critical radius from the CDM is also shown (dotted lines). Note that with increasing interface thickness the difference between the PFT and CDM predictions for the height of the nucleation barrier increases (see Table 1).
model potentials. The accuracy, however, will depend on how realistic the applied intermolecular potential is. Work is under- way to determine the thermodynamic and interfacial proper- ties of CO2 hydrate from MD simulations with realistic intermolecular potentials.
IV.2 Magnetic resonance imaging
The MRI intensity profile of water along the length of the par- tially water saturated sandstone core plug is shown in Fig. 10 as a function of time. Liquid CO2has been injected into the
Fig. 8 Snapshot of the water–CO2hydrate interface in molecular dynamics simulation (top). Note the regular clathrate cages inside the solid (on the right) and the distorted cages at the interface (on the left). For comparison, the ideal crystal structure (bottom left), and the real crystal with thermal fluctuations (bottom right) from MD simulations are also shown. We used SPC water and EPM2 CO2potentials.39The diameter of the H2O cages of circular projection (tetrakaidecahedra) is 0.866 nm. The rods at the center of the cages denote the CO2molecules. [These illustrations have been made using the Visual Molecular Dynamics package.41]
Table 1 Free energy of critical fluctuations (W*) vs.the interface thickness in the CO2hydrate system, as predicted by the phase field theory (PFT) and the classical droplet model (CDM) at 274 K and c¼0.033
d(nm) 1019WPFT/J 1019WCDM/J
0.5 3.18
1.0 2.66 3.76
1.5 2.23
Fig. 9 Height of the nucleation barrier for CO2hydrate formation as a function of interface thickness in the phase field theory (squares). For comparison, the prediction of the classical droplet model that assumes a sharp interface (d¼0) is also presented (circle). Note that ford!0 the PFT results converge to the value predicted by the CDM (dashed line: second order polynomial fitted to the PFT results).
plug at the face being at the left-hand side of this plot. The water saturation distribution, after 0.5 PV of liquid CO2dis- placed the water is shown by the upper saturation profile. At static conditions, the temperature was then lowered to 2C.
The downward arrow shows the trend of intensity variation with time. Decreasing intensity indicates the formation of hydrate since the MRI parameters were adjusted to image only the liquid hydrate precursors. Hydrate formation was indi- cated by a measurable decrease in the amount of free water and the simultaneous decrease in permeability with time. The noticeable difference in the change of the MRI intensity along the length of the plug suggests that the hydrate formation was primarily limited by the amount of CO2 available, i.e. the greatest change was near the CO2inlet face and the smallest of change was at the far end of the plug. The permeability decreased to zero as hydrate formed, implying that selected, but not readily detected, areas were completely filled with hydrate before all of the pore volume was filled.
A sequence of two-dimensional (2D) longitudinal images of a 4 mm thick slice along the center of the core plug during hydrate formation is shown in Fig. 11. These images imply that the hydrate primarily formed (areas of reduced intensity) near the CO2 inlet end of the plug, the right-hand side of these images. The hydrate formation process continued over a rela- tively long period of time, and the growth from right to the left is consistent with the initial nucleation of hydrate followed by the growth of hydrate crystals. These 2D images show that hydrate formed as a uniform band across the plug, rather than as separate, discrete clumps, again consistent with nucleation- initiated hydrate growth from the right-hand side of the plug.
V. Summary
We reported on advances made in two research areas related to the formation/dissolution of CO2 hydrate in aqueous solutions:
A theoretical development by us, based on the phase field theory of Gra´na´syet al.,16,17addresses the nucleation of CO2
hydrate in aqueous solution. First, we demonstrated for the hard-sphere system that this model is able to predict the nucleation rate quantitatively without adjustable parameters.
In contrast, the classical droplet model, used widely, predicts nucleation rates that are orders of magnitude too high. There- fore, we adapted the phase field theory to describe homoge- neous nucleation of hydrate from dissolved CO2 in aqueous solution. It has been demonstrated under typical conditions that the thickness of the hydrate-solution interface is compar- able with the size of nuclei, implying that the classical droplet model might be inaccurate. Indeed, as found for other systems, the nucleation barriers predicted by the phase field theory and
the classical droplet model differ considerably. Apparently, advanced models are needed to evaluate the rate of hydrate nucleation accurately. To improve the accuracy of predictions, work has been started to estimate the interface free energy and the width of the interface from molecular dynamics simula- tions for this system.
On the experimental side, we demonstrated the feasibility of magnetic resonance imaging as a tool for monitoring phase transitions involving hydrate in real porous media.
Acknowledgements
We thank R. L. Davidchack and B. B. Laird for communicat- ing us their new high-resolution density data for the (111), Fig. 11 A sequence of two-dimensional images that shows the forma- tion of CO2hydrate as function of time. The vertical axis displays posi- tion along the core height in cm. The horizontal axis is the length of the core plug. Top: when cooled to 2C. Center: After 30 min at 2C.
Lower: After 8 h at 2C. Gray scale indicates fraction of water not bound in hydrate.
Fig. 10 Formation of CO2 hydrate along the core as function of time.
(110), and (100) interfaces in the hard-sphere system prior to publication. This work has been supported by the Norwegian Research Council under project Nos. 153213/432 and 151400/
210, by the ESA under Prodex Contract No. 14613/00/NL/
SFe, and by the Hungarian Academy of Sciences under con- tract No. OTKA-T-037323. T. P. acknowledges support by the Bolyai Ja´nos Scholarship of the Hungarian Academy of Sciences.
References
1 L. Milich,Global Environ. Change, 1999,9, 179.
2 J. Mienert, K. Andreassen, J. Posewang and D. Lukas,Ann. N . Acad. Sci., 2000,912, 200.
3 W. J. Boettinger, J. A. Warren, C. Beckermann and A. Karma, Ann. Rev. Mater. Res., 2002,32, 163.
4 R. Kobayashi,Physica D, 1993,63, 410.
5 A. Karma and W.-J. Rappel,Phys. Rev. E, 1998,57, 4323.
6 J. A. Warren and W. J. Boettinger,Acta Metall. Mater., 1995, 43, 689.
7 W. L. George and J. A. Warren,J. Comput. Phys., 2002,177, 64.
8 A. Karma,Phys. Rev. E, 1994,49, 2245.
9 K. R. Elder, F. Drolet, J. M. Kosterlitz and M. Grant,Phys. Rev.
Lett., 1994,72, 677.
10 Nestler and A. A. Wheeler,Physica D, 2000,138, 114.
11 T. S. Lo, A. Karma and M. Plapp,Phys. Rev. E, 2001,63, 31 504.
12 F. Drolet, K. R. Elder, M. Grant and J. M. Kosterlitz,Phys. Rev.
E, 2000,61, 6705.
13 A. Graue, E. Aspenes, B. Kvamme, B. A. Baldwin, J. Stevens, D. P. Tobola and D. R. Zornes,Monitoring the Formation and Dissociation of Gas Hydrate in Reservoir Rock Using Magnetic Resonance Imaging, to be presented at the 2003 International Symposium of the Society of Core Analysts, Pau, France, Sept.
21–24, 2003.
14 R. L. Davidchack and B. B. Laird,J. Chem. Phys., 1998,108, 9452.
15 D. W. Oxtoby,Annu. Rev. Mater. Res., 2002,32, 39.
16 L. Gra´na´sy, T. Bo¨rzso¨nyi and T. Pusztai,Phys. Rev. Lett., 2002, 88, 206 105.
17 L. Gra´na´sy, T. Bo¨rzso¨nyi and T. Pusztai,J. Cryst. Growth., 2002, 237–239, 1813.
18 G. Caginalp and J. Jones,Ann. Phys. (NY), 1995,237, 66.
19 Y. C. Shen and D. W. Oxtoby, J. Chem. Phys., 1996, 105, 6517.
20 S. L. Wang, R. F. Sekerka, A. A. Wheeler, B. T. Murray, S. R. Coriell, R. J. Braun and G. B. McFadden, Physica D, 1993,69, 189.
21 J. W. Cahn and J. E. Hilliard,J. Chem. Phys, 1958,28, 258; 1959, 31, 688.
22 See for example K. F. Kelton,Solid State Phys., 1991,45, 75.
23 K. F. Kelton and A. L. Greer,Phys. Rev. B, 1988,38, 10 089.
24 K. R. Hall,J. Chem. Phys., 1972,57, 2252.
25 R. L. Davidchack and B. B. Laird,J. Chem. Phys., 1998,108, 9452.
26 We received new, high-resolution density data for the (111), (110), and (100) interfaces from R. L. Davidchack and B. B. Laird. They are fully consistent with their previous results for the (100) and (111) interfaces, published in ref. 25.
27 R. L. Davidchack and B. B. Laird,Phys. Rev. Lett., 2000,85, 4751.
28 B. Kvamme and H. Tanaka,J. Phys. Chem., 1995,99, 7114.
29 B. Kvamme, unpublished.
30 P. B. Stewart and P. Munjal,J. Chem. Eng. Data, 1970,15, 67.
31 D. Kashchiev and A. Firoozabadi, J. Cryst. Growth, 2002, 243, 476.
32 S. C. Hardy,Philos. Mag., 1977,35, 471.
33 R. M. Pratt, D.-H. Mei, T.-M. Guo and E. D. Sloan,J. Chem.
Phys., 1997,106, 4187.
34 H. Teng and A. Yamasaki,J. Chem. Eng. Data, 1998,43, 2.
35 L. Gra´na´sy and T. Pusztai,J. Chem. Phys., 2002,117, 10 121.
36 W. H. Shih, Z. Q. Wang, X. C. Zeng and D. Stroud,Phys. Rev. A, 1987,35, 2611.
37 P. R. ten Wolde, M. J. Ruiz-Montero and D. Frenkel,Phys. Rev.
Lett., 1996,75, 2714.
38 S. Auer and D. Frenkel,Nature, 2001,409, 1020.
39 J. G. Harris and K. H. Yung, J. Phys. Chem., 1995, 99, 12 021.
40 J. J. Hoyt, M. Asta and A. Karma,Phys. Rev. Lett., 2001,86, 5530.
41 W. Humphrey, A. Dalke and K. Schulten,J. Mol. Graphics, 1996, 14, 33.
