Solubility of carbon dioxide in water: some useful results for hydrate nucleation
Jesús Algaba Iván M. Zerón José Manuel Míguez Joanna Grabowska Samuel Blazquez Eduardo Sanz Carlos Vega Felipe J. Blas


AIP/123-QED

Solubility of carbon dioxide in water]Solubility of carbon dioxide in water: some useful results for hydrate nucleation



Laboratorio de Simulación Molecular y Química Computacional, CIQSO-Centro de Investigación en Química Sostenible and Departamento de Ciencias Integradas, Universidad de Huelva, 21006 Huelva Spain


Laboratorio de Simulación Molecular y Química Computacional, CIQSO-Centro de Investigación en Química Sostenible and Departamento de Ciencias Integradas, Universidad de Huelva, 21006 Huelva Spain


Laboratorio de Simulación Molecular y Química Computacional, CIQSO-Centro de Investigación en Química Sostenible and Departamento de Ciencias Integradas, Universidad de Huelva, 21006 Huelva Spain


Department of Physical Chemistry, Faculty of Chemistry and BioTechMed Center, Gdansk University of Technology, ul. Narutowicza 11/12, 80-233 Gdansk, Poland
Dpto. Química Física I, Fac. Ciencias Químicas, Universidad Complutense de Madrid, 28040 Madrid, Spain


Dpto. Química Física I, Fac. Ciencias Químicas, Universidad Complutense de Madrid, 28040 Madrid, Spain


Dpto. Química Física I, Fac. Ciencias Químicas, Universidad Complutense de Madrid, 28040 Madrid, Spain


Dpto. Química Física I, Fac. Ciencias Químicas, Universidad Complutense de Madrid, 28040 Madrid, Spain


Laboratorio de Simulación Molecular y Química Computacional, CIQSO-Centro de Investigación en Química Sostenible and Departamento de Ciencias Integradas, Universidad de Huelva, 21006 Huelva Spain
felipe@uhu.es



§ ABSTRACT 
In this paper, the solubility of carbon dioxide (CO_2) in water along the isobar of 400 bar is determined by computer simulations using the well-known TIP4P/Ice force field for water and TraPPE model for CO_2. In particular, the solubility of CO_2 in water when in contact with the CO_2 liquid phase, and the solubility of CO_2 in water when in contact with the hydrate have been determined. The solubility of CO_2 in a liquid-liquid system decreases as temperature increases. The solubility of CO_2 in a hydrate-liquid system increases with temperature. The two curves intersect at a certain temperature that determines the dissociation temperature of the hydrate at 400 bar (T_3). We compare the predictions with the T_3 obtained using the direct coexistence technique in a previous work. The results of both methods agree and we suggest 290(2) K as the value of T_3 for this system using the same cutoff distance for dispersive interactions. We also propose a novel and alternative route to evaluate the change in chemical potential for the formation of hydrate along the isobar. The new approach is based on the use of the solubility curve of CO_2 when the aqueous solution is in contact with the hydrate phase. It considers rigorously the non-ideality of the aqueous solution of CO_2, providing reliable values for driving force for nucleation of hydrates in good agreement with other thermodynamic routes used. It is shown that the driving force for hydrate nucleation at 400 bar is larger for the methane hydrate than for the carbon dioxide hydrate when compared at the same supercooling. We have also analyzed and discussed the effect of the cutoff distance of the dispersive interactions and the occupancy of CO_2 on the driving force for nucleation of the hydrate.


[
    Felipe J. Blas
    July 2024
==================














§ INTRODUCTION


At ambient conditions of temperature and pressure (298 K and 1 bar), the thermodynamically stable phase of water is the liquid phase. If the temperature is decreased at a constant pressure of 1 bar, the liquid is no longer the most stable phase and a first-order phase transition takes place at 273.15 K. Consequently, and according to the Thermodynamics laws, water must freeze. The new thermodynamically stable phase is the well-known ordinary ice, also known as Ih or hexagonal ice. This solid phase is formed by a crystalline structure characterized by the oxygen atoms forming hexagonal symmetry with nearly tetrahedral bonding angles. The same happens if the pressure is above ambient conditions up to 2100 bar, approximately. Above this pressure, water can freeze into other ices, including ice III, V, VI, among others, as the pressure is increased. <cit.> These are only some of the solid crystalline phases of the well-known polymorphic phases of water. However, this only happens if the original liquid phase is formed from pure water. When liquid water is mixed with another substance the story can be different.

There exist aqueous solutions of small compounds that exhibit different behavior when cooled down at constant pressure. Particularly, aqueous solutions of methane (CH_4), carbon dioxide (CO_2), nitrogen (N_2), hydrogen (H_2) or larger organic molecules, among many other different compounds, do not transform into a crystalline ice phase when the temperature is lowered. In fact, all these aqueous solutions freeze into new crystalline solid compounds named clathrate hydrates or simply hydrates. <cit.> Hydrates are non-stoichiometric crystalline inclusion compounds consisting of a network of hydrogen-bonding water molecules forming cages in which solutes (for instance, CH_4, CO_2, N_2 or H_2) are enclathrated at appropriate thermodynamic conditions of temperature and pressure.

Fundamental and applied research on hydrates and clathrates has been motivated by several reasons. First at all, hydrates are potential alternative sources of energy since huge amounts of CH_4 have been identified in hydrate deposits, either in the sea floor or in the permafrost frozen substrates, but their exploitation is not technically accessible yet due to a poor physicochemical characterization and various engineering issues. <cit.> Another remarkably relevant aspect of hydrates from both the scientific point of view and practical interest is the possibility to capture <cit.> and store CO_2. <cit.> This places gas hydrates at the center of environmental concerns regarding atmospheric greenhouse gases. Sequestration and capture of CO_2 in hydrates constitute a technological breakthrough which is seen as a promising alternative to other conventional methodologies for CO_2 capture, such as reactive absorption using amines and selective adsorption using adsorbent porous materials including sieves and zeolites. <cit.>

 

It is clear from the previous discussion that an accurate knowledge of  the thermodynamics and kinetics of the formation and growth of hydrates is necessary from the fundamental and practical points of view.  The thermodynamics of hydrates has been relatively well-established experimentally for years. <cit.> In addition, it is also possible to describe theoretically the phase equilibria of hydrates using the van der Waals and Platteeuw (vdW&P) formalism. <cit.> This approach, combined with an equation of State (EOS) allows us to satisfactorily determine the phase equilibrium of both pure hydrates and mixtures. <cit.> Additionally, from the point of view of molecular simulation, there has been an enormous development in techniques and methodologies for the study of the formation and dissociation of a huge variety of hydrates. <cit.> Particularly, several research groups have determined the phase equilibrium of CO_2<cit.> and CH_4 hydrates under oceanic crust conditions<cit.> using the direct coexistence technique. The precise knowledge of phase equilibria of hydrates, and particularly their phase boundaries, is essential to provide a detailed description of the kinetic and nucleation processes of these systems.

Unfortunately, a complete description from a molecular perspective of the mechanisms of growth and hydrate formation is far from being satisfactory. In the last few years, some of the authors of this work have been working on the development and use of the Seeding Technique, <cit.> in combination with the Classical Nucleation Theory (CNT), <cit.> to deal with several systems including the hard-sphere and Lennard-Jones models and more complex systems such as water and salty water. <cit.> More recently, we have extended the study to deal with methane hydrates. <cit.> It is important to recall here that Molinero and collaborators used the Seeding Technique to estimate nucleation rates of hydrates <cit.> modeled through the well-known mW water model. <cit.> Other authors have also contributed significantly to the understanding of the dynamics of nucleation and dissociation of hydrates from computer simulation. <cit.> This work constitutes the extension of our most recent study <cit.> on methane hydrates to deal with CO_2 hydrates. Before undertaking nucleation studies of CO_2 hydrates it is necessary to account for several issues, including the solubility of CO_2 in the aqueous solution when it is in contact with the CO_2-rich liquid phase and with the hydrate, an accurate prediction of the dissociation temperature, and the driving force for the nucleation.

The phase behavior of the CO_2 + water binary mixture is dominated by a large region of liquid-liquid (L_w-L_CO_2) immiscibility.^51,52 Since the critical point of pure CO_2 and a liquid-liquid-vapor (L_w-L_CO_2-V) three-phase line are at conditions similar to those at which the CO_2 hydrates are found, between 270-295 K and 10-5000 bar, approximately, another three-phase coexistence line involving a hydrate phase (i.e., a triple point that occurs at a certain temperature T_3 for each pressure) exhibits two branches. This is contrary to what happens with methane hydrates, that only exhibit one branch. <cit.> Fig. <ref> shows the pressure-temperature (PT) projection of the phase diagram of the CO_2 + water binary mixture. At pressures below 44.99 bar, the dissociation line is a H-L_w-V three-phase line at which the hydrate, the aqueous solution of CO_2, and the vapor phases coexist. Above that pressure, the hydrate and the solution coexist with a CO_2 liquid phase and a three-phase H-L_w-L_CO_2 line where the hydrate, the aqueous solution of CO_2, and the liquid phase of CO_2 coexist starts. Both branches meet at a Q_2 quadruple point located at 283 K and 44.99 bar (black filled circle) at which the hydrate, the aqueous solution, the CO_2 liquid, and the vapor phases coexist, <cit.> as can be seen in Fig. <ref>. Note that at Q_2 the L_w-L_CO_2-V three-phase line also meets with another H-L_CO_2-V three-phase in which the hydrate, the CO_2 liquid, and the vapor phases coexist at lower temperatures. In addition to this, there exists another quadruple point Q_1, located at 273 K and 12.56 bar (black filled square), at which the hydrate, the Ih ice, the solution, and the vapor phases coexist. This quadruple point connects the H-L_w-V three-phase line with a new three-phase H-Ih-V line involving the hydrate, the Ih ice, and the vapor phases that runs towards lower temperatures and pressures.




In this work, we concentrate on 400 bar of pressure (see the 400 bar isobar in Fig. 1 along which all the simulations are performed). At these conditions, the key solubility curves in the context of nucleation of CO_2 hydrates are the solubility of CO_2 in water when the aqueous solution is in contact with the CO_2 liquid phase and with the hydrate phase. In the first case, the solubility of CO_2 increases as the temperature is decreased. In the second case, as it occurs for the methane hydrate, there is little or no information from computer simulations or experiments. Here we determine the solubility of CO_2 in water from the hydrate along the isobar of 400 bar. This will allow us to estimate the dissociation line of the hydrate at this pressure, as we have done in our previous work for the case of the methane hydrate. <cit.>

The dissociation line of the CO_2 hydrate has been already determined by us several years ago. <cit.> It is important to mention that also other authors have obtained similar results using computer simulations <cit.> and free energy calculations. <cit.> Our previous results are slightly different than those found by Costandy and coworkers and Waage an collaborators since unlike dispersive interactions between water and CO_2 are different. However, we follow our previous work <cit.> and determine the dissociation line of the hydrate using the solubility curve of CO_2 in the aqueous solution when it is in contact with the CO_2 liquid phase and the hydrate. We have found that the new estimations agree with the initial prediction of Míguez et al. <cit.> within the corresponding uncertainties.

The formation of the CO_2 hydrate can be viewed as a chemical reaction in which water and CO_2 molecules “react” in the aqueous solution phase to form hydrate molecules. <cit.> The change in chemical potential of this reaction is the driving force for nucleation, Δμ_N. It is difficult to get good estimates of Δμ_N from experiments since it requires accurate values for a number of thermodynamic properties, including enthalpies and volumes of reactions, among other magnitudes. <cit.> Here we use the three independent routes introduced in our previous paper <cit.> to deal with the nucleation driving force for the nucleation of CO_2 hydrates. Particularly, we calculate the driving force for nucleation with respect to the state on the H-L_w-L_CO_2 three-phase line at 400 bar. Note that this point is well above the two quadruple points Q_1 and Q_2 shown in Fig. <ref>. In addition to this, we also propose a novel and alternative thermodynamic route based on the use of the solubility curve of CO_2 with the hydrate. This new route, that considers rigorously the non-ideality of the aqueous solution of CO_2 and provides reliable results of the driving force for nucleation, can be also used to determine Δμ_N of other hydrates.

The organization of this paper is as follows: In Sec. II, we describe the methodology used in this work. The results obtained, as well as their discussion, are described in Sec. III. Finally, conclusions are presented in Sec. IV.






§ METHODOLOGY


We use the GROMACS simulation package <cit.> to perform MD simulations. Computer simulations have been performed using three different versions of the NPT or isothermal-isobaric ensemble. For pure systems which exhibit fluid phases (pure water and pure CO_2) and aqueous solutions of CO_2 which exhibit bulk phases, we use the standard isotropic NPT ensemble, i.e., the three sides of the simulation box are changed proportionally to keep the pressure constant. For the hydrate phase, we use the anisotropic NPT ensemble in which each side of the simulation box is allowed to fluctuate independently to keep the pressure constant. This ensures that the equilibrated solid phase has no stress and that the thermodynamic properties are correctly estimated. The same ensemble is used to simulate the two-phase equilibrium between the hydrate and the aqueous solution of CO_2 (SL coexistence). Finally, the two-phase equilibrium between the solution and the CO_2 liquid phase is obtained using the NP_z𝒜T ensemble in which only the side of the simulation box perpendicular to the LL planar interface is allowed to change, with the interface area kept constant, to keep the pressure constant. For simulations involving LL and SL interfaces, we have used sufficiently large values of interfacial areas 𝒜. The thermodynamics and interfacial properties obtained from simulations of LL interfaces do not show a dependence on the surface area for systems with 𝒜>10× 10 σ^2. <cit.> Here σ is the largest Lennard-Jones diameter of the intermolecular potentials. In all simulations, the 𝒜 values used are higher than this value for LL and SL interfaces. See Sections III.A and III.C for the particular values used in this work.


In all simulations, we use the Verlet leapfrog<cit.> algorithm with a time steps of 2 fs. We use a Nosé-Hoover thermostat, <cit.> with a coupling time of 2 ps to keep the temperature constant. In addition to this, we also use the Parrinello-Rahman barostat <cit.> with a time constant equal to 2 ps to keep the pressure constant. We use two different cutoff distances for the dispersive and coulombic interactions, r_c=1.0 and 1.9 nm.
We use periodic boundary conditions in all three dimensions. The water-water, CO_2-CO_2, and water-CO_2 long-range interactions due to coulombic forces are determined using the three-dimensional Ewald technique. <cit.> Particularly, the real part of the coulombic potential is truncated at the same cutoff as the dispersive interactions. The Fourier term of the Ewald sums is evaluated using the particle mesh Ewald (PME) method. The width of the mesh is 0.1 nm, with a relative tolerance of 10^-5. In some calculations, we also use the standard long-range corrections for the LJ part of the potential to energy and pressure with r_c=1.0 nm. Water molecules are modeled using the TIP4P/Ice model <cit.> and the CO_2 molecules are described using the TraPPE model. <cit.> The H_2O–CO_2 unlike dispersive energy value is given by the modified Berthelot combining rule, ϵ_12=ξ(ϵ_11 ϵ_22)^1/2, with ξ=1.13. This is the same used by Míguez et al, <cit.> which allows us to predict accurately the three-phase hydrate–water–carbon dioxide coexistence or dissociation line of the CO_2 hydrate, particularly the coexistence temperature at the pressure considered in this work, 400 bar (see Fig. 10 and Table II of the work of Míguez and co-workers for further details). Very recently, we have demonstrated that the same molecular parameters are able to predict accurately the CO_2 hydrate-water interfacial free energy. <cit.>

Finally, uncertainties are estimated using standard deviation of the mean values or sub-block average method. Particularly, bulk densities in the LL and SL coexistence studies are obtained by averaging the corresponding density profiles over the appropriate regions sufficiently away from the interfacial regions. The statistical uncertainties of these values are estimated from the standard deviation of the mean values. Solubilities of CO_2 in all the liquid phases are calculated as molar fractions from the densities of both components and the corresponding errors are obtained from propagation of uncertainty formulae. Uncertainties associated to LL interfacial tension values, molar enthalpies, and partial molar enthalpies are estimated using the standard sub-block average method. Particularly, the production periods are divided into 10 (independent) blocks. The statistical errors are estimated from the standard deviation of the average.





§ RESULTS





 §.§ Solubility of carbon dioxide in water from the CO_2 liquid phase


We first concentrate on the solubility of CO_2 in the aqueous solution when the system exhibits LL immiscibility. In this case, there exists a coexistence between the water-rich and CO_2 liquid phases. We have used the direct coexistence technique to determine the solubility of CO_2 in the aqueous solution from the CO_2 liquid phase at several temperatures along the 400 bar isobar. Particularly, we have performed MD NP_z𝒜T simulations to ensure that temperature and pressure are constant. According to this, the planar interfacial area 𝒜=L_x× L_y is kept constant and only L_z is varied along each simulation. Here, L_x, L_y, and L_z are the dimensions of the simulation box along the x-, y-, and z-axis, respectively. In this work, the z-axis is chosen to be perpendicular to the planar interface. The initial simulation box is prepared in the following way. We build a slab of 2800 water molecules in contact, via a planar interface, with a second slab of 1223 molecules of CO_2. The dimensions of L_x and L_y of all the simulation boxes used in this part of the work are kept constant with L_x=L_y=3.8 nm (𝒜≃ 12× 12 σ^2). Since the pressure is constant, L_z varies along each simulation for all the temperatures considered. In this work, L_z varies from 11.06 to 12.29 nm. Simulations to calculate solubilities are run during 100 ns. The first 20 ns are used to equilibrate the system and the last 80 ns are used as the production period to obtain the properties of interest. We have also determined the LL interfacial tension and details of the simulations are explained later in this section.

Figure <ref> shows the density profiles of water and CO_2 as obtained by MD NP_z𝒜T simulations at 400 bar and temperatures from 250 to 310 K. For a better visualization for the reader, we only plot half of the profiles corresponding to one of the interfaces exhibited by the system. The right side of the figure corresponds to the CO_2 liquid phase and the left side to the aqueous liquid phase. We divide the inhomogeneous simulation box into 200 parallel slabs along the z–direction, perpendicular to the planar LL interface, to study the density profiles. Following the standard approach, density profiles are obtained assigning the position of each interacting site to the corresponding slab and constructing the molecular density from mass balance considerations.





As can be seen, density profiles of water (dashed curves) exhibit preferential adsorption at the interface at all temperatures. Particularly, water molecules are accumulated at the aqueous phase side of the interface. The relative maximum, which is identified with the accumulation of the water molecules, increases as the temperature of the system is decreased. The bulk density of water in the aqueous solution of CO_2 (left side of the figure) slightly decreases as the temperature is lower, especially in the range of 250-290 K. As can be seen, density profiles of water are nearly equal to zero at the bulk CO_2 liquid phase, indicating that solubility of water in that phase is completely negligible. Míguez et al. <cit.> have previously studied the LL interface of aqueous solutions of CO_2 at similar temperatures (287 and 298 K) but at lower pressures (P≤ 55 bar). These authors have found similar behavior for the density profiles of water but with an important exception: they exhibit the traditional shape of the hyperbolic tangent function in which water density decreases monotonically from the bulk density of water in the aqueous phase to zero in the CO_2 liquid phase.

The behavior and structure of the density profiles of CO_2 (continuous curves) along the interface are similar to those exhibited by other mixtures but with an important exception: CO_2 molecules exhibit activity on both sides of the liquid–liquid interface of the system. Particularly, there is an accumulation of CO_2 molecules at the CO_2 liquid phase side of the interface. This accumulation increases as the temperature of the system is decreased, as it happens with water molecules on the other side of the interface. The bulk density of CO_2 in both phases increases as the temperature is increased. This variation is more important in the CO_2 liquid phase (right side of Fig. <ref>). Contrary to what happens with water density in the CO_2 phase, the density of CO_2 in the aqueous solution is not negligible. This indicates that although the solubility of CO_2 in water is small (molar fraction of CO_2 between 0.04 and 0.09 in the range 310–250 K, respectively), its value is not so low as in the case of the solubility of water in CO_2. 

It is interesting to mention that density profiles of CO_2 also exhibit depressions in the aqueous solution side of the interface indicating desorption of CO_2 molecules in this region. The desorption of CO_2 molecules at the interface is correlated with the preferential adsorption of water molecules since the relative maxima and minima occur at the same position (z≈ -0.45 nm). Note that preferential adsorption and desorption of CO_2 molecules at the LL interface of this kind of aqueous solutions has not been previously seen in the literature. Particularly, Míguez et al. <cit.> only observe density profiles that exhibit preferential adsorption at the interface.




The solubilities of CO_2 in the aqueous phase have been determined from the information of the density profiles presented in Fig. <ref> at the corresponding temperatures.
Fig. <ref> shows the solubilities of CO_2 along the isobar at the temperatures considered. We have also included in the figure (inset) the same results obtained previously by us corresponding to the methane + water system. <cit.> As can be seen, the solubility decreases as temperature increases. This result is in agreement with our previous work in which we considered the solubility of methane in water along the same isobar (400 bar) and in contact with the gas phase (see the inset). <cit.> In this work, the study is firstly done using a relatively short cutoff distance for dispersive interactions, r_c=1.0 nm, which corresponds to a reduced cutoff value of r^*=r_c/σ≊ 3.16, with σ=0.31668 nm. Here σ is the length scale of the Lennard-Jones intermolecular interactions associated with the water model (TIP4P/ice). <cit.> In order to evaluate the effect of the cutoff distance, we have also determined the solubility of CO_2 using a larger cutoff value for the dispersive interactions (1.9 nm instead of 1.0 nm). As can be seen, the effect of increasing the cutoff is important in the whole range of temperatures considered. In particular, solubility increases between 17 %, at high temperatures (310 K), and 13 % at low temperatures (250 K). This is an expected result according to previous studies of the effect of the cutoff distance on fluid-fluid coexistence. <cit.>

We have checked that there is no a priori temperature limit to perform the simulations as the temperature is decreased. From this point of view, the solubility of CO_2 in the aqueous solution can be computed without any difficulty since we do not observe nucleation of the hydrate at low temperatures. This is in agreement with previous results obtained by Grabowska et al. <cit.> However, as the temperature is decreased the dynamics of the system slows down and the equilibration of the LL interface becomes more difficult and longer simulation runs are required to achieve equilibrated density profiles.


We have also determined the solubility using the standard long-range corrections to energy and pressure
to the Lennard-Jones part of the potential (dispersive interactions). According to our results,
although long-range corrections are able to improve the solubility results, differences between these results and those obtained with a cutoff of  1.9 nm  are still noticeable. In particular, the solubilities predicted using this approach are underestimated between 4 and 6 % along the isobar
at the temperatures simulated. It is interesting to compare the behavior of the CO_2 solubility, as a function of the temperature along the 400 bar, with that corresponding to methane obtained by us previously. <cit.> As can be seen in the inset of Fig. <ref>, the effect of the long-range correction on the dispersive interactions is slightly larger in the case of CO_2 than in methane. This is an expected result since the CO_2 molecules are modeled using three Lennard-Jones interaction sites and methane only with one, and also because the solubility of CO_2 is about ten times higher than the CH_4 solubility at the same thermodynamic conditions. It is also remarkable that the use of longer cutoff distances has a contrary effect on the solubility of CO_2 than in that of methane, i.e., the solubility of CO_2 increases with an increase of the cutoff whereas it decreases in the case of the solubility of methane. This is probably due to the presence of the quadrupolar moment of the CO_2 molecule and to the water–CO_2 interactions.

In the previous paragraphs, we have presented and discussed the results corresponding to the solubility of CO_2 in the aqueous solution when it is in contact with the CO_2 liquid phase. Since both phases are in contact through a planar LL interface and are in equilibrium at the same P and T, the chemical potential of water and CO_2 in both phases must satisfied that,


    μ_CO_2^I(P,T,x_CO_2^I)= μ_CO_2^II(P,T,x_CO_2^II)


and


    μ_H_2O^I(P,T,x_CO_2^I)= μ_H_2O^II(P,T,x_CO_2^II)


Here the superscripts I and II label the aqueous and the
CO_2 liquid phases, respectively. Note that we have expressed the chemical potential of water in each phase in terms of the corresponding CO_2 molar fractions. It is important to note that this is consistent from the thermodynamic point of view and it is always possible since we are dealing with a binary system that exhibits two-phase equilibrium. Following the Gibbs phase rule, such a system has two degrees of freedom, which according to Eqs. (<ref>) and (<ref>) are P and T. Consequently, the thermodynamic behavior of the system is fully described solving the previous equations since the composition of water in both phases, x_H_2O^I and x_H_2O^II, can be readily obtained as x_H_2O^I=1-x_CO_2^I and 
x_H_2O^II=1-x_CO_2^II.

According to our previous results shown in Fig. <ref>, the density of water in the CO_2 liquid phase is ρ_H_2O^II≈ 0, and consequently, x_H_2O^II=ρ_H_2O^II/(ρ_H_2O^II+ρ_CO_2^II)≈ 0
and x_CO_2^II=1-X_H_2O^II≈ 1.

Following the approximations of the previous paragraph combined with Eq. (<ref>), the chemical potential of CO_2 in the aqueous solution can be obtained from the chemical potential of pure CO_2 at the same P and T,


    μ_CO_2^I(P,T,x_CO_2^I)≈μ_CO_2^II(P,T,x_CO_2^II≈ 1)≈μ_CO_2^II(P,T)






The chemical potential of CO_2 along the isobar can be obtained from the thermodynamic relation,


    (∂(μ_CO_2/T)∂ T)_P,N_H_2O,N_CO_2=-h_CO_2T^2


where h_CO_2=h_CO_2(P,T) is the partial molar enthalpy of CO_2, and the derivative is performed at constant pressure, P, and number of water and CO_2 molecules, N_H_2O and N_CO_2, respectively. Since in our case the CO_2 liquid phase is essentially a pure CO_2 liquid, h_CO_2 is simply the molar enthalpy. Consequently, the chemical potential of CO_2, as a function of the temperature, along the 400 bar isobar can be obtained by integrating the Eq. (<ref>) as,


    μ_CO_2(T)k_BT=μ_CO_2(T_0)k_BT_0-_T_0^Th_CO_2(T')k_BT'^2 dT'


where k_B is the Boltzmann constant and T_0 is a certain reference temperature. Following our previous work <cit.>, we set μ_CO_2(T_0)=0. According to Eq. (<ref>), μ_CO_2(T) can be obtained by performing MD NPT simulations of pure CO_2 along the 400 bar isobar. In this case, since we are simulating a bulk phase, the standard NPT is used in such a way that the three dimensions of the simulation box are allowed to fluctuate isotropically. We use a cubic simulation box with 1000 CO_2 molecules. The dimensions of the simulation box, L_x, L_y, and L_z vary depending on the temperature from 3.9 to 4.2 nm. Simulations to calculate the molar enthalpy, at each temperature, are run during 100 ns, 20 ns to equilibrate the system and 80 ns as the production period to obtain h_CO_2. As in our previous work, we have not included the kinetic energy contribution (i.e. 5/2k_BT in the case of a rigid diatomic molecule such as CO_2). Note that this contribution is canceled out since we are evaluating chemical potential differences at constant P and T. In this work we choose as the reference temperature T_0=290 K. The reason for this selection will be clear later in the manuscript. Figure <ref> shows the chemical potential of CO_2 as a function of the temperature (blue curve). We have also included in the same figure the chemical potential values of the bulk methane taken from our previous work <cit.> (green curve) in order to compare both chemical potentials. Note that the reference temperature at which μ of the bulk methane is set to zero is 295 K.



 §.§ LL interfacial free energy


From the same simulations we have also obtained the LL interfacial tension, γ, from the
diagonal components of the pressure tensor. The vapour pressure corresponds to the normal component,
P≡ P_zz, of the pressure tensor. The interfacial tension is obtained using the well-known combination of the normal component and the tangential components, P_xx and P_yy through the mechanical route as, <cit.>


    γ=L_z/2[<P_zz>-<P_xx>+<P_yy>/2]


In Eq. (<ref>), the factor 1/2 reflects that during the simulations there exist two LL interfaces in the system, being L_z the size of the simulation box in the z direction perpendicular to the planar interface. Fig. <ref> shows the LL interfacial tension value as obtained from MD NP_z𝒜T simulations. Results obtained in our previous work corresponding to the LL interfacial tension of the methane + water mixture are also shown in the inset of the figure. <cit.> The interfacial tension decreases as the temperature is increased. We first calculate the interfacial tension using a cutoff value of 1.0 nm (green circles). In this case, we have used 50 ns for the equilibration period and 50 ns more for the production period in which the averages are calculated. Our results indicate that the values exhibit large fluctuations, especially at low temperatures. To improve our results, we have extended the simulations. The first 150 ns correspond to the equilibration period and the extra 150 ns are used to obtain the corresponding average values (production period). As can be seen (green diamonds), although the mean values obtained in both cases are similar, the error bars decrease, especially at the lowest temperature.



It is well-known that the equilibrium interfacial tension value associated with an interface critically depends on the molecular details. In particular, its value is very sensitive to the cutoff due to the dispersive interactions used during the simulations. <cit.> To account for the long-range interactions associated with the dispersive interactions, we have performed simulations using a cutoff distance of 1.9 nm, with 50 ns for equilibration and another 50 ns for production time (red circles). Note that in this case we are not using long-range corrections to energy and pressure. This value corresponds to a reduced cutoff distance r_c^*=r_c/σ=6, which is nearly the double value used in the first set of simulations. As can be seen, the main effect of increasing the cutoff distance is to decrease the interfacial tension values. Particularly, the effect is larger at low temperatures where the difference is about 3-4 mJ/m^2. However, at high temperatures, differences are about 1 mJ/m^2.

In order to account for the effect of the simulation length, we have extended the simulation at 310 K. As in the other set of simulations, we have equilibrated the system during the first 150 and used the next 150 ns to perform the corresponding averages (red triangles up). As can be seen, the new results are in practice identical to those obtained using only 50 ns for production time.

To be consistent with the calculations of the solubility of CO_2 in the aqueous solution, we have also determined the interfacial tension using the traditional LRC to energy and pressure. As in the previous case, we have only simulated using these corrections at 250 and 310 K (orange circles). As can be seen, this approximation is not able to provide consistent results at low temperatures when compared with the data obtained using r_c=1.9 nm. In particular, the interfacial tension is overestimated by more than 4.5 mJ/m^2, which represents more than 13% with respect to the value obtained using the larger cutoff distance. At the highest temperature, however, the overestimation of the interfacial tension is only about 0.5 mJ/m^2 (2%). Discrepancies at low temperatures between the results obtained using the largest cutoff distance (1.9 nm) and those using the traditional energy and pressure long-range corrections are probably due to the differences in densities in both liquid phases at these conditions.

We have also compared the predictions obtained from MD simulations with experimental data taken from the literature (magenta squares). <cit.> Unfortunately, to the best of our knowledge, there are no experimental data below 310 K. Our results are in good agreement with experimental measurements, although we slightly overestimate it by about 6%.


Finally, it is also interesting to compare the LL interfacial tension values obtained in this work for the CO_2 + water system with those obtained in our previous work <cit.> for the methane + water binary mixture. Computer simulation values of the former system are shown in the inset of Fig. <ref>. As can be seen, LL interfacial tension values of the methane + water system are approximately twice than those corresponding to the mixture containing CO_2. But perhaps the most interesting feature is that, as it happens with the solubilities (see Fig. <ref>), an increase of the cutoff distance due to the dispersive interactions has the opposite effect in the systems which containing methane and CO_2: an increase of the cutoff distance in the CO_2 system lowers the LL interfacial tension values, while in the methane mixture the interfacial tension values increase when the cutoff is larger (blue triangles up correspond to r_c=0.9 nm and dark green triangles down to r_c=1.7 nm). Also note that the effect of the long-range dispersive contributions is more important in the CO_2 + water mixture than in the system containing methane. As we have discussed previously, this could be due to the presence of the electric quadrupole of CO_2. 




 §.§ Solubility of carbon dioxide in  water from the hydrate phase


We have also determined the solubility of CO_2 in water when the aqueous solution is in contact with the hydrate along the 400 bar isobar at several temperatures. We first prepare a simulation box of CO_2 hydrate replicating a unit cell of hydrate four times along each spatial direction (4× 4 × 4), using 2944 water and 512 CO_2 molecules. This corresponds to a hydrate with the cages (8 cages per unit cell) fully occupied by CO_2 molecules. We equilibrate the simulation box for 40 ns using an anisotropic barostat along the three axes. This allows  the dimensions of the simulation box to change independently. The pressure is the same along the three directions and equal to 400 bar to allow the solid to relax and avoid any stress. In order to help the system to reach the equilibrium, we also prepare boxes of aqueous solutions with different concentrations of CO_2 depending on the temperature. This allows to reach the equilibrium as fast as possible in the last stage of the simulations when the hydrate and liquid phases are put in contact (see below). Particularly, the hydrate phase will grow or melt depending on the initial conditions, releasing/absorbing water and CO_2 molecules to/from the aqueous solution, until the solution phase reaches the equilibrium condition. Although the final state is independent of the initial CO_2 concentration in the aqueous phase, care must be taken in finite systems as those studied in this work. Initial conditions must be close enough to coexistence so that the system is able to reach equilibrium before exhaustion of any of the phases at coexistence. In this particular work, we have checked that density profiles of water and CO_2 in the aqueous phase reach the equilibrium value. This is practically done monitoring the averages profiles every 100 ns until no significant variations in their bulk region are observed. Once densities of water and CO_2 are obtained, the molar fraction of CO_2 in the aqueous solution is calculated from the corresponding averaged density values.

We use simulation boxes of solutions containing 4000 water molecules and varying the number of CO_2 molecules depending on temperature: 50 (250, 260, and 270 K), 120 (280 and 290 K), and 240 (295 K) CO_2 molecules. We equilibrate each simulation box during 40 ns using the isothermic-isobaric or NP_z𝒜T ensemble. In this case, two of the dimensions of the simulation boxes, arbitrarily named L_x and L_y, are kept constant (L_x=L_y vary between 4.77 and 4.82 nm (𝒜≃ 15× 15 σ^2) depending on the temperature) and equal to the values of two lengths of the simulation box of the hydrate. L_z is however allowed to vary to achieve the equilibrium pressure of 400 bar. Particularly, L_z varies from 10.28 to 10.53 nm depending on the temperature. Finally, the equilibrated hydrate and aqueous solution simulation boxes are assembled along the z direction sharing a planar solid-liquid interface with interfacial area 𝒜=L_x× L_y. We then perform simulations in the NPT ensemble using an anisotropic barostat with pressures identical in the three directions and equal to 400 bar. This allows the solid to relax and avoid any stress and obtain the correct value of the solubility at each temperature. Systems are equilibrated during 100 ns. After this, we run additional 300 ns to obtain the equilibrium density profiles of the system from 250 up to 295 K.


Fig. <ref> shows the density profiles of water and CO_2 molecules as obtained from anisotropic NPT simulations at 400 bar and temperatures ranging from 250 to 295 K. The density profiles have been obtained as explained in section III.A. Note that at temperatures above 295 K it is not possible to determine the solubility because the hydrate melts. In other words, there is a kinetic limit at high temperatures to determine the solubility of CO_2 from the hydrate.


As in the case of the LL coexistence described in section III.A, we only plot half of the profiles corresponding to one of the interfaces exhibited by the system. The right side of the figure corresponds to the hydrate phase and the left side to the aqueous solution phase. The density profiles in the hydrate phase exhibit the usual solid-like behavior for water and CO_2 molecules, with peaks at the corresponding crystallographic equilibrium position at which molecules are located in the hydrate. As can be seen, the density profiles at the lowest temperatures, from 250 up to 280 K show nearly the same structure, and only small differences are  observed at the hydrate-solution interface, as it is expected.

It is also interesting to analyze the behavior of profiles of water and CO_2 in the aqueous phase. The density profiles near the interface show some structural order due to the presence of the hydrate phase. Note that the positional order of the molecules is more pronounced at low temperatures, below T≤ 280 K. The bulk density of water (left side of the figure) slightly decreases as the temperature is increased, especially close to temperatures at which the hydrate melts. It is interesting to mention that bulk density profiles vary with temperature in the opposite way that when the aqueous solution is in contact with the CO_2 liquid phase (see Fig. <ref>).




The bulk density of CO_2 (in the aqueous solution phase) increases as the temperature is increased. From the inspection of Fig. <ref> it is clearly seen that the hydrate phase becomes less stable as the temperature approaches to 290-295 K. At lower temperatures, the hydrate-solution interface is located at z≈ 5 nm, approximately, with the hydrate phase showing 6-7 well-defined CO_2 layers. However, at 290 and 295 K only 5-6 layers can be observed in the hydrate phase, with the interface located at z≈ 6 nm. In fact, as we have previously mentioned, it is not possible to keep stable the hydrate at temperatures above 295 K, which eventually melts at higher temperatures.


We have calculated, from the information obtained from the density profiles, the solubility of CO_2 in the aqueous solution when it is in contact with the hydrate. As can be seen in Fig. <ref>, the solubility of CO_2 increases as the temperature is raised. We have also included in Fig. <ref> (inset) the same results obtained previously by us corresponding to the methane + water system. <cit.> It can be seen that our results are in agreement with the results of the solubility of methane. Contrary to what happens in the case of the methane hydrate, it is only possible to calculate the solubility up to a temperature of 295 K.
This is just a few degrees (around five) above the values of the three point temperature T_3 of the CO_2 hydrate at this pressure (see section III.D). In the case of our previous work, the hydrate was kept in metastable equilibrium at about 35 K over the dissociation temperature of the methane hydrate (295 K). <cit.>

We have also considered the impact of using different cutoff distances on solubilities. As in the case of the results presented in Section III.A, we have used three different cutoff distances for the dispersive interactions. Particularly, we use the same two values, 1.0 and 1.9 nm. We have also performed simulations using the standard long-range corrections to energy and pressure with a cutoff value of 1.0 nm. Contrary to what happens with the solubility of CO_2 in water when the solution is in contact with the other liquid phase (CO_2), the solubility does not depend on the cutoff distance, as can be seen in Fig. <ref>. Our results indicate that the long-range corrections due to the dispersive interactions have little or negligible effect on solubilities in aqueous solutions in contact with the hydrate. However, according to Figs. <ref> and <ref>, long-range interactions play a key role in the thermodynamic and interfacial properties in systems involving fluid phases. <cit.>

We have considered the effect of the CO_2 occupancy in the hydrate on the solubility in the aqueous solution. We have prepared the initial simulation boxes in a similar way to the case of full occupancy but with CO_2 occupying half of the small or D cages. Particularly, we use 2944 water (46× 4× 4 4) and 448 CO_2 molecules. This means that the occupancy of the large or T cages is 100% (384 CO_2 molecules) and the occupancy of the small or D cages is 50% (64 CO_2 molecules). This represents a 87.5% of occupancy of D and T cages. According to the experimental data, <cit.> the equilibrium occupancy of large or T cages of the CO_2 hydrate is nearly 100%. However, although there is a large discrepancy in measurements and predictions of the small cage occupancy, it is generally accepted that occupancy of small or D cages is approximately 30-60% depending on thermodynamic conditions. Notice that with this occupancy (i.e 87.5%) the ratio of water to CO_2 molecules in the hydrate is not 5.75 (as when the occupancy is 100%) but its value is now 46/7 ≈ 6.57. We follow the same procedure explained in the previous paragraph, with a cutoff distance of 1.0 nm. As a result, we obtain similar density profiles to those shown in Fig. <ref>. The solubility of CO_2 in water when the aqueous solution is in contact with the hydrate with 87.5% of occupancy is also shown in Fig. <ref>. As can be seen, the solubility of CO_2 when it is in contact with the hydrate with an occupancy of 87.5% is the same as that when it is in contact with the hydrate fully occupied within the error bars.

Finally, it is important to remark that, contrary to what we have found in our previous work for the solubility of methane in water, <cit.> we do not find a melting of the hydrate in a two-step process, i.e., a bubble of pure methane appears in the liquid phase as a first step and then the methane of the aqueous solution moves to the bubble and the methane from the hydrate moves to the aqueous solution as a second step. We find here that the hydrogen bonds of the layer of the hydrate in contact with the aqueous solutions break and the hydrate starts to melt. Particularly, when the temperature is increased the concentration of CO_2 in the aqueous solution increases in order to stabilize the hydrate phase. In the case of the methane hydrate, the amount of methane molecules releases to the aqueous phase to achieve the new equilibrium state is small (the solubility of methane in water is very small) and the metastable hydrate phase can exists above the T_3. However, in the case of the CO_2 hydrate, the CO_2 saturates the aqueous phase (the solubility of CO_2 in water increases greatly with the temperature). The hydrate becomes unstable, the hydrogen bonds of the hydrate layer next to the aqueous solution breaks, and the hydrate finally melts.







 §.§ Three-phases coexistence from solubility calculations


We have obtained the solubility of CO_2 in the aqueous solution, as a function of the temperature at a fixed pressure of 400 bar, when it is in contact with the CO_2 liquid phase (Section III.A) and with the hydrate (Section III.C). In both cases, the system exhibits two-phase coexistence. It is interesting to represent both solubilities in the same plot as we did in our previous study, <cit.> and as it is shown in Fig. <ref>. Since one of the solubility curves is a decreasing function of the temperature and the other an increasing function of the temperature, there exists a certain temperature, that we will call T_3 for reasons that will be clear soon, at which both solubilities are equal at 400 bar.

The points of the solubility curve of CO_2 in the aqueous solution from the CO_2 liquid phase correspond to thermodynamic states at which the pressure and the chemical potentials of water and CO_2 in the aqueous phase are equal to those in the CO_2 liquid phase. In addition to this, the points of the solubility of CO_2 in water from the hydrate phase correspond to states at which the pressure is also the same and at which the chemical potentials of both components in the aqueous phase are equal to those in the hydrate phase. Consequently, at T_3, the temperature, pressure, and chemical potentials of water and CO_2 in the aqueous solution, CO_2 liquid, and hydrate phases are the same. This means that the point at which the two solubility curves cross represents a three-phase coexistence state of the system at 400 bar. This is also known as the dissociation temperature of the CO_2 hydrate at the corresponding pressure (400 bar).

The value obtained in this work for the T_3 is 290(2) K when the occupancy of the hydrate is 100% (all the T and D cages are occupied by CO_2 molecules). We assume here an uncertainty of 2 K for the dissociation temperature of the hydrate, following the same estimation of the T_3 error of the methane hydrate determined in our previous work. <cit.> We have also determined the dissociation temperature of the hydrate when the occupancy of the small or D cages is 50% (87.5% overall occupancy) using r_c=1.0 nm. In this case, T_3=290.5(2) K, which is the same value obtained for the fully occupied hydrate within the error. Both dissociation temperature results seem to be occupancy-independent with the employed methodology and are in good agreement (within the corresponding uncertainties) with the value obtained by Míguez et al. using the direct coexistence technique,<cit.> 287(2) K. It is important to remark here that we are using the same models for water (TIP4P/ice)<cit.> and CO_2 (TraPPE),<cit.> the same unlike dispersive interactions between both components, and cutoff distance for dispersive interactions (r_c=1.0 nm) than in the work of  Míguez et al. <cit.> At this point it is important to remark that the system sizes of this work are different than that used in the work of Míguez et al., <cit.> and this could have a subtle effect in the T_3 because of finite-size effects, as it has been found for the melting point of ice Ih. <cit.> The experimental value of T_3 at 400 bar is 286 K so that the force field used in  this work provides a quite reasonable prediction. 

Other authors have determined the T_3  for this system from computer simulation. Costandy et al.<cit.> have calculated the dissociation temperature of the CO_2 hydrate at 400 bar using the direct coexistence technique. They obtained a value of 283.5(1.7) K. Although they also used the same water and CO_2 models, a number of differences lead to a slightly different value of T_3: different unlike dispersive interactions between water and CO_2 and cutoff distance for dispersive interaction (1.1 nm). Waage and collaborators <cit.> have also determined the dissociation line of the hydrate using free energy calculations. They also use the same models for water and CO_2 but different unlike dispersive interactions between them. In this case, the cutoff distance for dispersive interactions is r_c=1.0 nm. These authors calculate the dissociation temperature of the hydrate at 200 and 500 bar. The values obtained are 283.9(1.7) and 284.8(0.9) K, respectively. Interpolating to 400 bar, the T_3 is 284.5 K, in good agreement with the results of Costandy et al.<cit.> Unfortunately, the result obtained here can not be compared with the predictions of Costandy et al.<cit.> and Waage and collaborators <cit.> since dispersive interactions are not the same as those used here.

Finally, it is important to focus on the effect of the cutoff distance used to evaluate the long-range dispersive interactions. The T_3 value of 290(2) K has been obtained using a cutoff distance of 1.0 nm. We have also analyzed the solubilities of CO_2 from the CO_2 liquid and the hydrate phases using a much larger cutoff distance (i.e., 1.9 nm). As we have previously shown, the solubility of CO_2 from the CO_2 liquid phase increases when the value of the cutoff is increased. On the other hand, the solubility of CO_2 in the hydrate phase is not affected by the use of larger cutoff values. Consequently, the combined effect of the increase of the cutoff distance of the dispersive interactions is an increase of the T_3 since it is the intersection of the two solubility curves shown in Fig. <ref>. Particularly, the dissociation temperature of the hydrate is now found at 292(2) K, 2 K above the T_3 observed with a cutoff distance of 1.0 nm, approximately. 

It is interesting to compare the effect of the cutoff due to the dispersive interaction in both CO_2 and methane hydrates. As can be seen, in the case of the methane hydrate, T_3 is shifted towards lower temperatures, by 2 K, when the cutoff is increased. In the current case (CO_2 hydrate), we observed the opposite effect, i.e., T_3 increases when the cutoff distance is increased. This is the same effect as observed for the solubility curve of CO_2 in the aqueous solution in contact with the CO_2 liquid phase. This effect, contrary to that observed in the methane hydrate, could be due to the electrostatic interactions of the quadrupole of CO_2 with other CO_2 molecules and also with water molecules, that it is not present in the case of the methane. We think this issue deserves a more detailed study but this is out of the scope of the current work.

We have determined the dissociation line of the CO_2 hydrate at 400 bar from the calculation of the solubility of CO_2 when the aqueous solution is in contact with the other two phases in equilibrium, the CO_2 liquid phase and the hydrate phase. Grabowska and collaborators <cit.> have already demonstrated that this route allows to determine T_3 of hydrates. This work confirms that this methodology is a good alternative to the direct coexistence method. Particularly, it shows a slightly better efficiency compared with the other technique (lesser simulation times are required) and provides consistent values of T_3. 






 §.§ Driving force for nucleation of hydrates


The dissociation line of the CO_2 hydrate separates its phase diagram in two parts in which two different two-phase coexistence regions exist. <cit.> At a certain pressure, for instance, 400 bar, at temperatures above T_3 the system exhibits LL immiscibility between an aqueous solution and a CO_2 liquid phase. Note that the solubility of water in CO_2 is very small and the CO_2 liquid phase can be considered pure CO_2 in practice. However, at temperatures below T_3 the system exhibits SL phase equilibrium between the hydrate and a fluid phase (water or CO_2 depending on the global composition of the system). This is consistent with the nature of the dissociation or three-phase line at which the hydrate, aqueous solution, and CO_2 liquid phases coexist.

The fluid phase in equilibrium with the hydrate below T_3 depends on the global composition of the system. Here we assume that the hydrate is fully occupied by CO_2 molecules, i.e., 8 CO_2 molecules for every 46 water molecules according to the stoichiometry of hydrates type sI. Let be N_H_2O and N_CO_2 the number of water and CO_2 molecules used in the fluid phases during the simulations, respectively. If the ratio N_H_2O/N_CO_2> 5.75, one should have hydrate–water phase separation (below T_3). However, if N_H_2O/N_CO_2< 5.75, one should have a hydrate – CO_2 phase system for T<T_3.

As described by Kashchiev and Firoozabadi <cit.> and Grabowska and collaborators, <cit.> the formation of a hydrate in the aqueous solution phase can be viewed as a chemical reaction that takes place at constant P and T,








    CO_2 (aq,x_CO_2) +
    5.75 H_2O (aq,x_CO_2) 
    → [CO_2(H_2O)_5.75]_H



Since we work at constant pressure in this work (P=400 bar), we drop the dependence of P in the rest of equations. Assuming that all the cages of the hydrate are filled, a unit cell of CO_2 hydrate is formed by 46 water molecules and 8 CO_2 molecules, i.e., 1 CO_2 molecule per 46/8=5.75 water molecules. According to this, Eq. (<ref>) considers the hydrate as a new compound formed from one molecule of CO_2 and 5.75 molecules of water. We can also associate to this compound one unique chemical potential for the hydrate at T, μ_H^H(T). Note that this chemical potential is simply the sum of the chemical potential of CO_2 in the solid plus 5.75 times the chemical potential of water in the solid, i.e.,


    μ^H_H(T)=μ_CO_2^H(T)+5.75 μ_H_2O^H(T)


According to the previous discussion, the compound [CO_2(H_2O)_5.75]_H is simply the “hydrate” and we call one “molecule” of the hydrate in the solid the molecule [CO_2(H_2O)_5.75].

Following Kashchiev and Firoozabadi, <cit.> we denote the driving force for nucleation of the hydrate formed from the aqueous solution with a concentration x_CO_2 at T as,


    Δμ_N(T,x_CO_2)   =μ^H_H(T) 
        -μ^aq_CO_2(T,x_CO_2)-5.75 μ^aq_H_2O(T,x_CO_2)


Note that Δμ_N in this paper is Δμ_nucleation in our previous paper. <cit.> Δμ_N also depends on pressure but since we are working at constant pressure (400 bar), we drop the pressure dependence from all the equations in this study. μ_H^H(T) has been previously defined in Eq. (<ref>) as the chemical potential of the “hydrate molecule” in the hydrate phase, and μ^aq_CO_2(T,x_CO_2) and  μ^aq_H_2O(T,x_CO_2) are the chemical potentials of CO_2 and water in the
aqueous solution, respectively, at T and molar fraction of CO_2, x_CO_2. Note that the composition of CO_2 in Eq. (<ref>) is, a priori, independent of the pressure and temperature selected. In other words, one could have different driving forces for nucleation, at a given P and T, changing the composition of the aqueous solution (for instance, in a supersaturated solution of CO_2). However, there exists a particular value of x_CO_2 which is of great interest from the experimental point of view. Experiments on the nucleation of hydrates are performed when the water phase is in contact with the CO_2 liquid phase through a planar interface. Since both phases are in equilibrium at P and T, the solubility of CO_2 in water (molar fraction of CO_2 in the aqueous solution) is fully determined since x^eq_CO_2≡ x_CO_2^eq(T). Following the notation of Grabowska and coworkers, <cit.> the driving force for nucleation at experimental conditions is given by,


    Δμ^EC_N(T)   =μ^H_H(T) -μ^aq_CO_2(T,x_CO_2^eq(T)) 
        -5.75 μ^aq_H_2O(T,x_CO_2^eq(T))


Note that Δμ^EC_N depends only on T (and on P but in this work we are working at the same P=400 bar). We provide here valuable information for this magnitude when the molecules are described using the TIP4P/Ice and TraPPE models for water and CO_2, respectively.

In the next sections we concentrate on the driving force for nucleation at experimental conditions, Δμ^EC_N(T), obtained using four different routes. In the first one (route 1), we use the definition of the driving force of nucleation given by Eq. (<ref>). In the second one (route 2), we use the solubility curves of CO_2 from the hydrate and the CO_2 liquid phase. In third one (route 3), we use the enthalpy of dissociation of the hydrate and assume that it does not change with temperature nor composition. In the fourth one (route 4), we propose a novel methodology based on the use of the solubility curve of CO_2 with the hydrate, valid not only for the CO_2 hydrate but also for other hydrates. This route can be used to determine Δμ_N at any arbitrary temperature and mixture composition and not only at experimental conditions. Finally, we discuss the results obtained using the different routes and compare the driving force for nucleation of the CO_2 hydrate with that of the methane hydrate previously obtained by us in a previous work. <cit.> 




  §.§.§ Route 1 for calculating Δμ^EC_N.


Route 1 was proposed and described in our previous work <cit.> and we summarize here only the main approximations and the final expression of the driving force for nucleation. To evaluate Δμ^EC_N(T,x^eq_CO_2) in Eq. (<ref>) we need to calculate the chemical potential of the “hydrate molecule” in the hydrate phase, and the chemical potentials of CO_2 and water in the aqueous phase  at a supercooled temperature T below T_3. The change in the hydrate chemical potential when the temperature passes from T_3 to T can be evaluated in a similar way as that for pure CO_2 from T_3 to T. In fact, this later change has been already calculated in Section III.A using Eq. (<ref>) and evaluated from the corresponding thermodynamic integration using computer simulations in the NPT ensemble.

The chemical potential of water in the aqueous phase at T can be estimated using the procedure of Grabowska and collaborators <cit.> applied to the case of the CO_2 hydrate according to Eqs. (20)-(26) of their paper. This is done in two steps. In the first step, the change in the chemical potential of the solution when its temperature passes from T_3 to T is approximated by that of pure water calculated from thermodynamic integration (see Eqs. (24) and (26) of the work of Grabowska and collaborators). The second step involves the change in the chemical potential of water in the solution when the composition of CO_2 changes from x^eq_CO_2(T_3) to x^eq_CO_2(T) (see Eq. (25) of our previous work).  The rigorous calculation of this contribution requires the knowledge of the activity coefficient of water in an aqueous solution with a given composition of water, γ^aq_H_2O(T,x^eq_H_2O), at T and T_3. Grabowska and coworkers assume that this magnitude, in the case of an aqueous solution of methane, γ^aq_H_2O≈ 1 since the solution is very diluted. We follow here the same assumption.

Using these approximations it is possible to compute the driving force at experimental conditions for nucleation given by Eq. (<ref>). The final expression is given by,





    Δμ^EC_N(T,x^eq_CO_2)k_BT=
    -_T_3^Th_H^H(T')- {h_CO_2^pure(T')+5.75
    h_H_2O^pure(T')}k_BT'^2dT'
    - [k_BTln{x^eq_H_2O(T)}-k_BT_3ln{x^eq_H_2O(T_3)}]




Here h_H^H=H/N_CO_2 is the enthalpy H of the hydrate per CO_2 molecule and N_CO_2 is the number of CO_2 molecules in the hydrate. Note that Eq. (<ref>) is consistent with the view of Kashchiev and Firoozabadi <cit.> of the hydrate as a new compound formed from one molecule of CO_2 and 5.75
 molecules of water when the hydrate is fully occupied. Also note that it is necessary to use that the driving force for nucleation at T_3 is equal to zero, i.e.,


    Δμ^EC_N(T_3,x^eq_CO_2)   =μ^H_H(T_3) -μ^aq_CO_2(T_3,x^eq_CO_2(T_3)) 
        -5.75 μ^aq_H_2O(T_3,x^eq_CO_2(T_3))=0


This is equivalent to arbitrary set to zero the chemical potentials of CO_2 and water in the hydrate at T_3. Eq. (<ref>) is similar to Eq. (<ref>) (route 3 or dissociation route) but taking into account two effects: (1) the temperature dependence of molar enthalpies of hydrate, CO_2, and water; and (2) the change of composition of the solution when passing from T_3 to T (see route 3 below for further details).


As we have previously mentioned, each change in the chemical potentials needed to compute the driving force for nucleation is obtained by evaluating molar enthalpies of pure CO_2, water, and hydrate phases involved in the integrals given by Eq. (<ref>). Note that the chemical potential of CO_2 in the aqueous solution has been already obtained in Section III.A. Particularly since in this route one is interested in computing Δμ_N at experimental conditions, the chemical potential of CO_2 in the aqueous solution is equal to that of pure CO_2 at the same P and T. Consequently, it has been obtained from the integration of the molar enthalpy at several temperatures according to Eq. (<ref>).


In the case of water, the chemical potential can be obtained performing MD NPT simulations of pure water along the 400 bar isobar. As in the case of pure CO_2, we also use the standard NPT is used in such a way that the three dimensions of the simulation box are allowed to fluctuate isotropically. We used a cubic simulation box with 1000 H_2O molecules. The dimensions of the simulation box, L_x, L_y, and L_z vary depending on the temperature from 3.14 to 3.08 nm. Simulations to calculate the molar enthalpy, at each temperature, are run during 100 ns, 20 ns to equilibrate the system, and 80 ns as the production period to obtain h^pure_H_2O.

In the case of the pure hydrate, we have obtained the chemical potential in a similar way, performing simulations in the NPT ensemble using an isotropic barostat at 400 bar. At the beginning of each simulation, we use a cubic box formed by 27 replicas of the unit cell in a 3× 3× 3 geometry. The dimensions of the simulation box vary between 3.85 and 3.62 nm depending on the temperature. As in the rest of simulations, we calculate the enthalpy at different temperatures, from 260 to 295 K. Simulations are run during 100 ns, 20 ns to equilibrate the system and 80 ns to calculate the molar enthalpy of the hydrate.




  §.§.§ Route 2 for calculating Δμ^EC_N.


Route 2 was also proposed and described in our previous work <cit.>. This route is inspired by the work of Molinero and coworkers<cit.> and we summarize here only the main approximations and the final expression of the driving force for nucleation. According to this, it is possible to find a different, but an equivalent, thermodynamic route to calculate the driving force for the nucleation of methane hydrates. We check in this work whether this approach can also be used to deal with CO_2 hydrates. Let us consider Eq. (<ref>) at experimental conditions, i.e., at the equilibrium composition of CO_2 in the aqueous solution when it is in contact via a planar interface with a CO_2 liquid phase (L), x^eq_CO_2(T|L). Note that the vertical line represents flat interface equilibrium with
the liquid CO_2. To clarify the derivation of the final expression, we write explicitly the solubility of water in the solution, at experimental conditions, as x^eq_H_2O(T|L). Obviously, as we have mentioned previously in Section III.A, the solubility of water in the solution can be obtained readily as  x^eq_H_2O(T|L)=1-x^eq_CO_2(T|L). 

We now assume that the chemical potential of the ideal solution's components can be expressed, in general, in terms of the chemical potentials of the pure components in the standard state and their molar fractions. In other words, since the molar fraction of CO_2 in the solution is small, we are assuming that water is the dominant component (solvent) and the CO_2 is the minor component (solute) in the mixture. Under these circumstances, the activity coefficients of water and CO_2 are close to one. <cit.> According to this and following our previous work,<cit.> Eq. (<ref>) can be written as,


    Δμ^EC_N(T)   =
    -k_BTln[x^eq_CO_2(T|L)x^eq_CO_2(T|H)] 
       -5.75k_BTln[x^eq_H_2O(T|L)x^eq_H_2O(T|H)]


x^eq_CO_2(T|H) and x^eq_H_2O(T|H) represent the molar fraction of CO_2 and H_2O in the solution when it is in equilibirum via a planar interface (vertical line) with the hydrate phase (H), respectively. Note that Δμ^EC_N(T) and all the molar fractions also depend on pressure but in this work we work at fixed P=400 bar. This is the equation obtained previously by us considering the driving force for the nucleation of the methane hydrate.<cit.> As we will see later in this section, Eq. (<ref>) does not provide reliable values for the driving force of nucleation of the CO_2 hydrate, contrary to what happens with the methane hydrate. The solubilities of methane in the solution when is in contact with the methane phase and with the hydrate are one order of magnitude lower  than those of CO_2 in the case of the CO_2 hydrate. Consequently, this route can be useful only in cases in which the solubility of the guest is extremely low.




  §.§.§ Route 3 (dissociation) for calculating Δμ^EC_N.


It is possible to estimate the driving force for nucleation of a hydrate using a simple and approximate route based on the knowledge of the enthalpy of dissociation of the hydrate. <cit.> 
The dissociation enthalpy of the hydrate, h_H^diss, is defined as the enthalpy change of the process, <cit.>


    [CO_2(H_2O)_5.75]_H→CO_2 (liq) +
    5.75 H_2O (liq)


Dissociation enthalpies are usually calculated assuming that the hydrate dissociates into pure water and pure CO_2. Note that this corresponds to the definition of enthalpy of dissociation and that in reality CO_2 will be dissolved in water and an even smaller amount of water will be dissolved in the CO_2 liquid phase. We have determined the dissociation enthalpy of the hydrate simply by performing simulations of the pure phases (hydrate, water, and CO_2) at several temperatures at 400 bar.

According to our previous work, <cit.> we evaluate the driving force for nucleation assuming the following approximations: (1) the enthalpy of dissociation of the hydrate, h_H^diss does not change with the temperature; (2) its value can be taken from its value at T_3; and (3) enthalpy of dissociation does not vary with composition of the aqueous solution containing CO_2 when the temperature is changed. According to this, Δμ_N^EC is given by,


    Δμ_N^EC=k_BT
    _T_3^Th^diss_Hk_BT'^2 dT'≊ -h^diss_H (T_3) (1-TT_3)


Note that Eq. (<ref>) reduces to Eq. (<ref>) under the approximations used in this route. 



  §.§.§ Route 4 for calculating Δμ^EC_N.



The driving force for nucleation of the CO_2 hydrate, at any arbitrary temperature, T_N, and molar fraction of CO_2 in the aqueous solution, x_CO_2^(N), at fixed pressure is defined as,


    Δμ_N(T_N,x^N_CO_2)   =μ^H_H(T_N) -μ^aq_CO_2(T_N,x^N_CO_2) 
        -5.75 μ^aq_H_2O(T_N,x^N_CO_2)


Note that Δμ_N also depends on pressure. However, since we work at constant pressure (P=400 bar), we drop the pressure dependence from equations from this point. It is also important to recall that since we are assuming that all cages of the hydrate are filled, the chemical potential of a “hydrate molecule” in the hydrate phase depends only on temperature. Finally, the chemical potentials of CO_2 and water also depend on the molar fraction of water in the aqueous solution, x_H_2O^N. Since we are dealing with a binary mixture, x_H_2O^N=1-x^N_CO_2. For simplicity, we choose x^N_CO_2 as independent variable of the chemical potentials of CO_2 and water in the solution. 




The driving force for nucleation of the CO_2 hydrate, Δμ_N, depends on T_N and x_CO_2^N and both are independent variables. This means that route 4 is valid for calculating the driving force for nucleation at any T_N and x_CO_2^N. As we will see later, the method can be particularized to evaluate Δμ_N at experimental conditions. In this case, Δμ_N=Δμ_N^EC(T_N)=Δμ_N^EC(T_N,x_CO_2^eq(T_N|L)), as we have previously mentioned.

To evaluate Δμ_N we need to calculate the chemical potential of the “hydrate molecule”, μ^H_H(T_N), at a supercooled temperature T_N, and the chemical potentials of CO_2 and water molecules of an aqueous solution of CO_2 with molar fraction x_CO_2^N at the same temperature. This route is based on the use of the solubility curve of the hydrate with temperature, at constant pressure, previously described in Section III.C. A schematic depiction of the curve and the thermodynamic route for obtaining Δμ_N at arbitrary T_N and x_CO_2^N is presented in Fig. <ref>. Let us consider a reference state in our calculations at temperature T_ref on the solubility curve in contact with the hydrate. As it will be clear later, the particular value of T_ref is not important since we are dealing with differences of chemical potentials and the final value of Δμ_N does not depend on the election of the reference state. Due to this, the reference state does not appear in Fig. (<ref>).

The first contribution to Δμ_N in Eq. (<ref>) is the chemical potential of the “hydrate molecule” in the hydrate phase at T_N. The chemical potential μ_H^H can be obtained using the Gibbs-Helmholtz thermodynamic relation for pure systems,


    (∂(μ_H^H/T)∂ T)_P, N_H=-h_H^HT^2


where h_H^H=h_H^H(T) is the molar enthalpy of the “hydrate molecule” and the derivative is performed at constant pressure, P, and number of “hydrate molecules”, N_H. Note that here h_H^H represents the enthalpy of the hydrate per molecule of CO_2 according to the definition in Section III.E.1. The chemical potential of the “hydrate molecule” in the hydrate phase at a supercooling temperature T_N can be obtained by integrating the Eq. (<ref>) from T_ref to T_N as,


    μ_H^H(T_N)k_BT_N=
    μ_H^H(T_ref)k_BT_ref
    -_T_ref^T_Nh_H^H(T)k_BT^2 dT


where k_B is the Boltzmann constant.

The last two contributions to the driving force for nucleation in Eq. (<ref>), μ^aq_CO_2 and μ^aq_H_2O, need to be evaluated at temperature T_N
and molar fraction x_CO_2^N. Individual chemical potentials of CO_2 and water in the solution at a given temperature are not easy to evaluate, as we have seen in routes 1 and 2. However, it is possible to use the solubility curve of CO_2 with the hydrate to overcome this problem. 
 
Let be T_i the temperature at which the aqueous solution with molar fraction x_CO_2^N is in equilibrium with the hydrate phase as indicated in Fig. <ref>. Since both phases are in equilibrium at these conditions, the chemical potentials of CO_2 and water in the hydrate phase and in the aqueous solution are equal,


    μ_CO_2^H(T_i)=
    μ^aq_CO_2(T_i,x_CO_2^N)



    μ_H_2O^H(T_i)=μ^aq_H_2O(T_i,x_CO_2^N)


Note that x_CO_2^N=x^eq_CO_2(T_i|H) according to the nomenclature used in Section III.E.2 (route 2) and in our previous paper. <cit.> The vertical line here represents that the aqueous solution is in equilibrium with
the solid hydrate via a flat interface.

 Combining Eqs. (<ref>) and (<ref>) with Eq. (<ref>), that gives the chemical potential of the “hydrate molecule” in terms of the chemical potentials of CO_2 and water in the hydrate phase, we obtain,


    μ_H^H(T_i)=μ_CO_2^aq(T_i,x_CO_2^N)
    +5.75 μ_H_2O^aq(T_i,x_CO_2^N)


Eq. (<ref>) is the heart of route 4. According to it, the combination μ_CO_2^aq+5.75 μ_H_2O^aq is known along the solubility curve of the hydrate at any temperature T_i: it is equal to the chemical potential of the "hydrate molecule" at the temperature considered. This apparently simple result allows to calculate accurately the driving force for nucleation at any temperature and composition of the solution using a one-step thermodynamic integration. As it will be clear at the end of this section, this method can be used to determine the driving force for nucleation of other hydrates. 

In the first step, we calculate the difference of μ_CO_2^aq
+5.75 μ_H_2O^aq between the reference state (ref) at T_ref and a second state (i) at T_i, both on the solubility curve of CO_2 with the hydrate as indicated in Fig <ref>. According to Eq. (<ref>), this is completely equivalent to evaluate the difference of μ_H^H between T_ref and T_i along the solublity curve. This change can be evaluated using again the Eq. (<ref>) (Gibbs-Helmholtz relation) and integrating between the two temperatures,


    μ_H^H(T_i)k_BT_i=
    μ_H^H(T_ref)k_BT_ref
    -_T_ref^T_ih_H^H(T)k_BT^2 dT



In the second step, that involves the difference between the chemical potentials of CO_2 and water in solution at temperatures T_i and T at constant molar fraction x_CO_2^N, Δμ_CO_2^aq and Δμ_H_2O^aq, can be obtained from the Gibbs-Helmholtz equation for CO_2 and water,


    (∂(μ_α^aq/T)∂ T)_P,x_α=-h_α^aqT^2







Here α={CO_2, H_2O} and represents one of the components of the mixture. Note that the partial derivative is calculated at constant composition. In this case, the composition corresponds to that of the aqueous solution in equilibrium with the hydrate phase at T_i. h_α^aq is the partial molar enthalpy of component α in the aqueous solution. The partial molar enthalpy is defined as,


    h_α^aq=N_A(∂ H∂ N_α)_P,T,N_βα=
    lim_Δ N_α→ 0N_A(Δ HΔ N_α)_P,T,N_βα


where N_A is the Avogadro's number and H is the aqueous solution's enthalpy. The limit can be numerically evaluated computing the enthalpy for two systems that have the same number of water molecules and different number of CO_2 to evaluate the partial molar enthalpy of CO_2. The partial molar enthalpy of water can be estimated in a similar way, i.e., the number of molecules of CO_2 in the system is kept constant while the number of water molecules changes. According to this, it is possible to evaluate the variation of the chemical potential of CO_2 and water from T_i to T_N, Δμ_CO_2^aq  and Δμ_H_2O^aq, from the knowledge of the partial molar enthalpies of both components. In particular, the combination of the chemical potentials of CO_2 and water, as a function of T_N, can be obtained by integrating Eq. (<ref>) as,




    μ_CO_2^aq(T_N,x_CO_2^N)+5.75
    μ_H_2O^aq(T_N,x_CO_2^N)k_BT_N=
    μ_CO_2^aq(T_i,x_CO_2^N)+5.75
    μ_H_2O^aq(T_i,x_CO_2^N)k_BT_i-_T_i^T_Nh_CO_2^aq(T,x_CO_2^N)+5.75
    h_H_2O^aq(T,x_CO_2^N)k_BT^2dT




Now, it is possible to find a closed expression for evaluating Δμ_N
at arbitrary T_N and x_CO_2^N in terms of the enthalpies of the "hydrate", CO_2, and water molecules. Using Eqs. (<ref>), (<ref>), (<ref>), and (<ref>), the driving force for nucleation can be written as,




    Δμ_N(T_N,x^N_CO_2)k_BT_N=
    -_T_i^T_Nh_H^H(T)- {h_CO_2^aq(T,x_CO_2^N)+5.75
    h_H_2O^aq(T,x_CO_2^N)}k_BT^2dT







We recall here that h_H^H is the enthalpy of the “hydrate molecule” per molecule of CO_2. Since we are assuming that the hydrate is fully occupied, the factor in the partial molar enthalpy of water must be 46/8=5.75 to be consistent with the stoichiometry of the unit cell. It is important to remark two important aspects of this expression. As we have previously mentioned, Δμ_N does not depend on the reference state. Note that the two integrations of the molar enthalpy of the “hydrate molecule” between T_ref and T_N, given by Eq. (<ref>), and between T_ref to T_i, given by Eq. (<ref>), are now expressed as a single integration of the molar enthalpy of the “hydrate molecule“ between T_i and T_N. In other words, since the driving force for nucleation does not depend on the reference state, the initial state of Eq. (<ref>) is simply T_i. This is also related with another important fact: the driving force for nucleation of the hydrate is zero not only at T_3 but also along the whole solubility curve of the hydrate.



The second interesting aspect of Eq. (<ref>) is that the integrand of the right term can be formally written as an enthalpy of dissociation of the hydrate that depends on T_N and x^N_CO_2. Under this perspective, Eq. (<ref>) resembles Eq. (<ref>) of route 3 since it has the same mathematical form.





the integrand of the right term can be rewritten in terms of the enthalpy of dissociation of the hydrate, as h_H^diss=h_H^H-(h_CO_2^aq+5.75 h_H_2O^aq). Note that h_H^diss depends on T_N and x_CO_2^N. A closed expression of the driving force for nucleation can be expressed in a more compact way as,



    Δμ_N(T_N,x^N_CO_2)k_BT_N=
    -_T_i^T_Nh_H^diss(T,x^N_CO_2)k_BT^2dT



Eq. (<ref>) resembles Eq. (<ref>) of route 3. In Eq. (<ref>), h_H^diss depends on T and x^N_CO_2 and Δμ_N is calculated assuming explicitly that the hydrate dissociates into water and CO_2 in the aqueous solution at these conditions. However, Eq. (<ref>) assumes that h_H^diss is constant and water and CO_2 dissociates into pure water and pure CO_2.







Eq. (<ref>) is a rigorous and exact expression (within the statistical uncertainties of the simulation results) obtained only from thermodynamic arguments for calculating the driving force for nucleation of the CO_2 hydrate at any T_N and x_CO_2^N. Obviously, this route is general and can be used to calculate driving forces for nucleation of other hydrates from the knowledge of the solubility curve of the corresponding guest with the hydrate.

Let us now apply this route to the particular case of the CO_2 hydrate and evaluate Δμ_N^EC(T) at experimental conditions, i.e., with x_CO_2^N≡ x_CO_2^eq(T|L). Each of the chemical potential changes in Eqs. (<ref>) or (<ref>) can be obtained evaluating the molar enthalpy of the “hydrate molecule” and the partial molar enthalpies of CO_2 and water in the aqueous solution in Eq. (<ref>). In the case of the pure hydrate, the change in the chemical potential can be obtained performing simulations in the NPT ensemble using an isotropic barostat at 400 bar. At the beginning of each simulation, we use a cubic box formed by 27 replicas of the unit cell in a 3× 3× 3 geometry. The dimensions of the simulation box vary between 3.85 and 3.62 nm depending on the temperature. As in the rest of simulations, we calculate the enthalpy at different temperatures, from 260 to 295 K. Simulations are run during 100 ns, 20 ns to equilibrate the system and 80 ns to calculate the molar enthalpy of the hydrate.




As can be seen in Fig. <ref>, the solubility curve of CO_2 with the liquid phase is a convex and decreasing function of temperature. According to this, it is not possible to reach this curve from the solubility curve of CO_2 with the hydrate from a temperature T_i below T_3 following the one-step thermodynamic path of route 4 (see also Fig. <ref>). A feasible solution is to choose a T_i value above T_3 where the hydrate–solution coexistence line is metastable. However, the T_3 of the CO_2 hydrate is located at T_3=290 K (with molar fraction x_CO_2^eq(T_3)=0.05), and the solubility curve with the hydrate can be only obtained up to 295 K, as it is shown in Fig. <ref>. This means that Eq. (<ref>) can be only applied to calculate Δμ^EC_N for supercoolings above 270-280 K, approximately (depending on the value of the cutoff). In the case of the methane hydrate it is possible to calculate the hydrate–solution equilibrium curve at temperatures significantly higher than T_3. This allows to evaluate Δμ^EC_N at lower tempertatures than in the case of the CO_2 hydrate (see the inset of Fig. <ref> in this work and Fig. 4 of our previous work <cit.>).

To overcome this problem, we propose to use Eq. (<ref>), for several values of T_i<T_3, and to perform an extrapolation of the composition at the temperature at which we evaluate Δμ^EC_N, as it is indicated schematically in Fig. <ref>. According to this, μ_CO_2^aq and μ_H_2O^aq, or the combination of both according to Eq. (<ref>), can be obtained by performing NPT MD simulations of the solution along the 400 bar isobar at constant composition. As in the case of the “hydrate molecule”, since we are simulating bulk phases, the standard NPT is used in such a way that the three dimensions of the simulation box are allowed to fluctuate isotropically. We evaluate the partial molar enthalpies of both components at five different concentrations: x_CO2=0.0680, 0.0521, 0.0413, 0.0335, and 0.0215. These values are the compositions of the aqueous solution along the solubility curve of CO_2 with the hydrate, x_CO_2^eq(T_i|H), at T_i=295, 290, 285, 280, and 270, respectively. We calculate numerically the derivative of Eq. (<ref>) by computing the enthalpy for two different systems that have the same number of water molecules and different number of CO_2 molecules and the same number of CO_2 molecules and different number of water molecules to determine h^aq_CO_2 and h^aq_H_2O, respectively. Particularly, h^aq_H_2O is obtained from the difference of the enthalpies of the aqueous solution using 990 and 1010 water molecules (Δ N_H_2O=20) for all the temperatures and compositions of the mixtures. In the case of h^aq_CO_2, we have used different number of CO_2 molecules depending on the composition of the mixture: 20 and 24 (Δ N_CO_2=4) for x_CO2=0.0215, 32 and 38 (Δ N_CO_2=6) for x_CO2=0.0335, 40 and 46 (Δ N_CO_2=6) for x_CO2=0.0413, 50 and 60 (Δ N_CO_2=10) for x_CO2=0.0521, and 68 and 78 (Δ N_CO_2=10) for x_CO2=0.0680. Simulations to calculate the enthalpy, at each temperature, are run during 300 ns, 50 ns to equilibrate the system and 250 ns as the production period. Dividing the enthalpy difference, Δ H, by the difference of the number of CO_2 molecules, Δ N_CO_2, and of the number of water molecules, Δ N_H_2O, and multiplying by Avogadro's constant, we get estimations of h^aq_CO_2 and h^aq_H_2O at several compositions of the mixture and temperatures.





Fig. <ref> shows the partial molar enthalpy of CO_2 at five constant compositions, as functions of the temperature. The composition of the mixture in each curve corresponds to the molar fraction of the solution at several temperatures, T_i, along the solubility curve of CO_2 with the hydrate, as indicated in the previous paragraph. The difference of the chemical potential of CO_2 along the integration path (see the Fig. <ref>), Δμ_CO_2^aq, is represented in the inset of Fig. <ref>. Note that we have five different curves, one corresponding to each temperature value T_i. We follow the same approach and represent the partial molar enthalpy of water at the same compositions in Fig. <ref>. In the inset is now depicted 5.75Δμ_H_2O^aq, as function of the temperature, which represents the change of the rest of the "hydrate molecule" dissociated in the solution according to the nomenclature used in this section. As can be seen, partial molar enthalpies of CO_2 and water show small variations with the composition and decrease as the temperature decreases. The variation of both chemical potentials also exhibit similar behavior, i.e., they decrease as the temperature is lowered when keeping the composition of the mixture constant. In addition, each Δμ_CO_2^aq→ 0 and Δμ_H_2O^aq→ 0 at the temperature associated with the corresponding composition on the solubility curve of CO_2 with the hydrate.



Using the values of μ_CO_2^aq and 5.75 μ_H_2O^aq, as functions of temperatures and at the molar fractions of CO_2 considered previously, it is possible to evaluate Δμ_N according to the scheme indicated in Fig. <ref>. Figure <ref> shows Δμ_N, as a function of x_CO2, at four temperatures below T_3 (250, 260, 270, and 280 K). As can be seen, Δμ_N follows a linear dependence with the composition of the solution. We have performed linear regressions using the five values of x_CO_2, for each temperature, to obtain Δμ^EC_N at experimental conditions. We have also represented the compositions of the solution in equilibrium with the CO_2 liquid phase on the solubility curve of CO_2, x^eq_CO_2(T|L) at the four temperatures (crosses and stars). As can be seen, we have two different values of 
x^eq_CO_2(T|L) depending on the cutoff due to the dispersive interactions used, r_c=1.0 (crosses) and  1.9 nm (stars).



Δμ_CO_2^aq+5.75 Δμ_H_2O^aq at a given temperature T and at the appropriate composition of the solution in equilibrium with the CO_2 liquid phase, x^eq_CO_2(T|L), can be obtained in a second step, as indicated in Fig. <ref>, according to the following approach. The values of Δμ_CO_2^aq and 5.75 Δμ_H_2O^aq, at any T<T_3, are obtained using Eqs. (<ref>) or (<ref>), with α=CO_2 and H_2O, and taking into account that μ_CO_2^aq+5.75 μ_H_2O^aq at each T_i is equal to μ_H^H(T_i) according to Eq. (<ref>).





It is also interesting to show μ_CO_2^aq+5.75 μ_H_2O^aq, as a function of composition, evaluated at several temperatures below the T_3 of the hydrate as shown in Fig. <ref> (see also Fig. <ref>). Note that we have set to zero the chemical potentials of CO_2 and water in the hydrate at T_3=290 K. This representation contains the same information as Fig. <ref> but allows to discuss important aspects related with the approximations of route 3. We have represented μ_CO_2^aq+5.75 μ_H_2O^aq at T=250, 260, 270, and 280 K at the compositions x^eq_CO_2 previously selected (symbols). In addition to that, we have also represented the value of the change in the chemical potential of the hydrate, at the four temperatures, in equilibrium with the aqueous solution along the solubility curve of the hydrate. Note that the value of μ_CO_2^aq+5.75 μ_H_2O^aq at 270 K is equal to that of the hydrate, -11.2 k_BT, since it is in equilibrium with the solution with x^eq_CO_2=0.0215. The same is true at 280 K, i.e., μ_CO_2^aq+5.75 μ_H_2O^aq=μ_H^H=-5.4 k_BT, state at which the solution with x_CO_2=0.0335 is in equilibrium with the hydrate.



 
The values of μ_CO_2^aq+5.75 μ_H_2O^aq, at the corresponding temperatures, follow a linear dependence with the composition of the solution. We have performed linear regressions using the five values obtained in this work, for each temperature, and the corresponding lines are also shown in Fig. <ref>. We have also represented the compositions of the solution in equilibrium with the CO_2 liquid phase on the solubility curve of CO_2, x^eq_CO_2(T|L) at the four temperatures (crosses and stars). Note that these values fit in an excellent way to the linear regression. As can be seen, we have two different values of 
x^eq_CO_2(T|L) depending on the cutoff due to the dispersive interactions used, r_c=1.0 (crosses) and  1.9 nm (stars). See Fig. <ref> and Section III.D for further details. According to this, it is possible to know accurately the values of μ_CO_2^aq+5.75 μ_H_2O^aq by extrapolating the linear fits depending on the temperature and the composition (see Fig. <ref>). These values, in combination with the values of μ_H^H(T) (already obtained in routes 1 and 2), can be used to predict with confidence the driving force for nucleation of the CO_2 hydrate using this new approach. Note that results obtained from Fig. <ref> are the same than those obtained from Fig. <ref> since they contain the same information.

Before finishing this section, Fig. <ref> contains valuable information that deserves to be discussed in detail. Dashed curves represent the values of μ_CO_2^aq+5.75 μ_H_2O^aq obtained at the lowest concentration considered, x_CO2=0.0215, and assuming that the variation with composition follows the approximation used in route 2 proposed by us in our previous work<cit.> and also used by Molinero and co-workers <cit.> and by us in this work. In other words, the difference of the chemical potentials of CO_2 and water when the composition is varied can be calculated assuming that the activity coefficients of both component are close to 1 or are similar in the solution. Unfortunately , this approximation is not valid for CO_2 hydrates. As can be seen in Fig. <ref>, this approach (route 2) underestimates the value of μ_CO_2^aq+5.75 μ_H_2O^aq more than 0.8 k_BT (using r_c=1.0 nm and more than 1 k_BT using r_c=1.9 nm) at 250 K with respect to the value obtained from route 4. As we will see later in the next section, this result explains from a thermodynamic perspective why the route 2 can not be used with confidence to estimate the driving force for nucleation of the CO_2 hydrate.




  §.§.§ Evaluation of Δμ^EC_N using different routes.




We have obtained the driving force for nucleation of the CO_2 hydrate using the four routes presented in the previous sections. All the results have been obtained using a cutoff distance for the dispersive interactions r_c=1.0 nm. As we have seen in the previous sections, this value of r_c gives a dissociation temperature of the hydrate T_3=290(2) K. The results obtained using the different routes are presented in Fig. <ref>. The route 1 given by Eq. (<ref>) predicts an almost linear behavior of Δμ_N^EC with the temperature. This results is in agreement with our previous results obtained for the driving force for nucleation of the methane hydrate. <cit.>

We have also used the novel route proposed in this work (route 4), based on the use of the solubility curve of CO_2 with the hydrate, given by Eq. (<ref>). As we have mentioned in the previous section, the route 4 should provide reliable values of Δμ_N^EC since it is based on rigorous thermodynamic integration calculations. The only approximation made are the extrapolations of Δμ_N to the x^eq_CO_2 on the solubility curve of CO_2 with the liquid at the corresponding temperatures. However, we think this is a good approach taken into account the low values of the concentration and the results presented in Figs. <ref> and <ref>. As can be seen in Fig. <ref>, small differences are seen between results obtained from routes 1 and 4. Route 1 slightly underestimates the driving force for nucleation in nearly the whole range of temperatures considered in this work, especially in the intermediate range of temperatures.





As in our previous work, <cit.> we have also used the dissociation route (route 3) according to Eq. (<ref>), proposed by Kashchiev and Firoozabadi. <cit.> The agreement between the results from route 3 and 4 is good, especially at low supercoolings. The dissociation route overestimates the values of Δμ_N^EC obtained from the route 4 about 0.1 k_BT at 260 K, approximately. This represents a value 6.5% higher than the value obtained from the route 4, the maximum difference found between both approaches in the whole range of temperatures.




Finally, we have also obtained Δμ_N^EC, as a function of the temperature, via the route 2 proposed by us in our previous work <cit.> and inspired by the work of Molinero and collaborators. <cit.> This route, given by Eq. (<ref>), entails crude approximations. As can be seen, the route 2 is not able to provide reliable predictions of Δμ_N^EC in the whole range of supercoolings. In fact, it overestimates its value by 0.32 k_BT, a value 20% higher than that obtained using the more rigorous route 4 in a wide range of temperatures. This result is in agreement with the findings observed in Figs. <ref> and <ref> and it is a direct consequence of the main approximation made in route 2: that the activity coefficients of water and CO_2 are equal to one. Although this is a good approximation for the methane hydrate, <cit.> it is not a realistic option for the CO_2 hydrate. The root of this behavior must be found in the large differences in solubility of methane in water compared with that of 
CO_2 in water (in contact with both, the gas/liquid phase and the hydrate phase). See the Figs. <ref> and <ref> and the corresponding insets.




  
In summary, the route 2 is not in general a good choice for calculating driving forces for nucleation of hydrates. This route can be used when the solubility of guest molecules in water is extremely low. The route 3 is an easy and fast way to estimate Δμ_N^EC values. However, we do not recommend this route in general except for cases in which solubility of the guest molecules in water is extremely low, as in the case of route 2. Finally, the route 4 proposed here in the most rigorous and nearly exact way to evaluate driving forces for nucleation of hydrates.


It is important to discuss the effect of the cutoff distance due to the dispersive interactions on Δμ_N^EC. We have already analyzed the effect of r_c on the solubility curve of CO_2 in contact with both, the CO_2 liquid (Fig. <ref>) and the hydrate (Fig. <ref>). Although the solubility curve of CO_2 with the hydrate is practically unaffected when r_c is changed from 1.0 to 1.9 nm, the situation is completely different for the solubility curve of CO_2 with the liquid. As a consequence of this, the T_3 of the CO_2 hydrate changes from 290(2) K when r_c=1.0 nm to 292(2) K when r_c=1.9 nm. Obviously, this change must also affect to the values of Δμ_N^EC. We have obtained the driving force for nucleation following the route 4 using a cutoff distance r_c=1.9 nm and results are compared with those using r_c=1.0 nm. As can be seen in Fig. <ref>, the main effect is to displace the curve towards higher temperatures. This is an expected result due to the difference in the T_3 values using different cutoff distances. However, it is clearly seen that differences between both curves increases as the temperature is decreased: the difference between both values in absolute value is 0.125 k_BT at 290 K, approximately, but that difference increases up to 0.286 k_BT at 260 K, approximately. This is more than double of the value of the difference predicted at 290 K, suggesting that the increase of the cutoff distance has a deep effect on the driving force for nucleation of the system.



To check the real impact of the cutoff distance of the dispersive interaction on the driving force for nucleation, we have plot Δμ_N^EC, as a function of the supercooling Δ T, instead of the absolute temperature T. This allows to compare both results at the same supercooling and to have a clearer picture of this effect. As can be seen in Fig. <ref>, there is an important effect on Δμ_N^EC when r_c is changed from 1.0 to 1.9 nm. For instance, at |Δ T|≈ 25 K, Δμ_N^EC changes from -1.5 to -1.65 k_BT when the cutoff distance is increased. According to this, the driving force for nucleation of the hydrate is 10% larger when r_c=1.9 nm than that obtained using 1.0 nm. This effect is not negligible. The origin of this displacement is due to the strong dependence of the solubility of CO_2 in water on r_c (aqueous solution in contact with the CO_2 liquid phase). Due to the effect of the cutoff distance on Δμ_N^EC, appropriate values of r_c are required in order to obtain reliable values of this magnitude.

It is also very interesting to compare the driving force for the nucleation of the CO_2 and methane hydrates at the same pressure. We have determined in our previous work <cit.> the driving force at experimental conditions along the same isobar. We also present these results in Fig. <ref>, obtained using the route 1 and a cutoff distance of r_c=0.9 nm. We compare these results with that obtained for the CO_2 hydrate using a cutoff distance of 1.0 nm. As can be seen, the driving force of the CO_2 hydrate is, in absolute value, lower than that of the methane hydrate along the isobar of 400 bar. For instance, at a supercooling of |Δ T|≈ 25 K, the driving force for the methane hydrate is Δμ^EC_N≈ -1.7 k_BT. At the same supercooling for the CO_2, Δμ^EC_N≈ -1.5 k_BT. This means that the driving force for the nucleation of the CO_2 hydrate, at 400 bar, is 13% lower than that of the methane hydrate at the same supercooling (Δ T=25 K). According to this, the nucleation of the methane hydrate should be more favorable than that of the CO_2 hydrate. Obviously, this would be true if the other factors that affect the nucleation rate of the hydrates are equal, i.e., the water–hydrate interfacial energy.

We have also considered the effect of the occupancy of CO_2 in the hydrate on the driving force for nucleation. Particularly, we study hydrates with 7 CO_2 molecules per unit cell, i.e., 50% of occupancy in the small or D cages and 100% of occupancy in the large or T cages, which is equivalent to 87.5% of overall occupancy. According to the work of Kashchiev and Firoozabadi, <cit.> the formation of a hydrate in the aqueous solution phase can be described as the chemical reaction of Eq. (<ref>). This reaction can be viewed as the formation of a “hydrate molecule” per each CO_2 molecule in the aqueous solution.

However, since we now calculate and compare driving forces for nucleation of hydrates with different occupancies, it is more convenient to write Eq. (<ref>) per cage of hydrate formed from the aqueous solution than per CO_2 molecule used to form the hydrate from the solution. In the case of a hydrate fully occupied, the reaction is the same in both descriptions since an unit cell of hydrate is formed from 8 cages (6 T and 2 D cages) and it is occupied by 8 CO_2 molecules as well. Let's define the occupancy, x_occ as the fraction of cages occupied by the CO_2, x_occ=n_CO_2/n_cg, where n_CO_2 and n_cg are the number of CO_2 molecules and cages per unit cell. When the occupancy is 100%, x_occ=8/8=1 and x_occ=7/8=0.875 when 
it is 87.5%. According to this, the formation of one cage of hydrate, with occupancy 87.5%, from the aqueous solution phase can be viewed as a classical chemical reaction that takes place at constant P and T,



    x_occ CO_2 (aq) +
    5.75 H_2O (aq) 
    ⟶ [(CO_2)_x_occ(H_2O)_5.75]_H







In this particular case, since each unit cell of CO_2 hydrate is formed from n_cg=8 cages and 46 water molecules, we only need 7/8=0.875 CO_2 molecules (i.e., an occupancy x_occ=7/8=0.875) and 46/8=5.75 water molecules in the solution to form 1 cage of hydrate with the desired occupancy (7 CO_2 molecules per unit cell). The compound [(CO_2)_x_occ(H_2O)_5.75]_H is simply a “cage” of hydrate. According to this, we call [(CO_2)_x_occ(H_2O)_5.75] a “molecule” of one cage of the hydrate in the solid. Note that stoichiometry of reaction given by Eq. (<ref>) is in agreement with an unit cell of this partially occupied hydrate, formed from 8 “cages” of hydrate with 8× 0.875=7 CO_2 molecules and 8× 5.75=46 water molecules.




We have used the route 1 described in Section III.E.1 with a cutoff distance for the dispersive interactions of r_c=1.0 nm. We have followed the same procedure previously explained but instead of simulating a hydrate fully occupied by CO_2 molecules we have considered a hydrate with occupancy of the small or D cages of 50% (87.5% overall occupancy of the hydrate). To be consistent with the description of the previous paragraph, we have used Eq. (<ref>) to evaluate the driving force for nucleation of the partially occupied hydrate per cage of hydrate instead of per CO_2 molecule. According to this, the corresponding molar enthalpy of the hydrate, h_H^H, as a function of the temperature, must be expressed as an enthalpy per cage of the hydrate, h̃_H^H,



    h_H^H=HN_CO_2=
    HN_cg(N_cg/N_CO_2)=
    h̃_H^H(N_cg/N_CO_2)=
    h̃_H^Hx_occ


Here H is the enthalpy of the hydrate, and N_CO_2 and N_cg are the total number of CO_2 molecules and cages used in the simulations, respectively.   
N_CO_2=n_cells× n_CO_2 and N_cg=n_cells× n_cg, with n_cells=3× 3× 3=27 the number of unit cells in the simulations. Note that the enthalpy per cage, h̃_H^H, is obtained dividing the enthalpy of the hydrate by the total number of cages, N_cg=n_cells× n_cg=27× 8=216. H is calculated in the same way as in Section III.E.1 for the fully occupied hydrate but now using a simulation box with 189 CO_2 molecules and 27 replicas of the unit cell in a 3× 3× 3 geometry used in Section III.E.1 for the 1242 water molecules.

Fig. <ref> shows the enthalpy per cage of the hydrate h̃_H^H partially occupied by CO_2 molecules as a function of temperature (blue diamonds). We also present the enthalpy per cage of the hydrate with 100% of occupancy (red circles). Note that 
h̃_H^H is equal to the enthalpy of the hydrate per molecule of CO_2, h_H^H, in the case of full occupancy. As can be seen, the enthalpy per cage when occupancy is 87.5% is systematically less negative than h̃_H^H when the hydrate is fully occupied. The difference between both values is ∼ 3 kJ/mol. This is an expected result since there is one CO_2 molecule less per unit cell (7 instead of 8) in the hydrate with occupancy of 87.5%. Although the difference is below 1%, there is less CO_2-water favorable (negative) dispersive interactions and this contributes to increase the energy and consequently the enthalpy of the system. Note that the lattice parameters of the unit cell (for a certain P and T) of the hydrate depend on the occupancy. Particularly, it becomes about 0.16% smaller when the occupancy changes from 100% to 87.5%.


Once h̃_H^H(T) is known, it is possible to use Eq. (<ref>) to evaluate the driving force for nucleation of the hydrate. However, Eq. (<ref>) is only valid for hydrates with 100% occupancy. It is possible to reformulate route 1 for hydrates partially occupied taken into account that the enthalpy of the hydrate is expressed per cage of the hydrate, h̃_H^H, and using the appropriate stiochiometry when the hydrate has an occupancy of x_occ. According to this, the driving force for nucleation per cage of hydrate is given by,





    Δμ^EC_N(T,x^eq_CO_2)k_BT=
    -_T_3^Th̃_H^H(T')- {x_occ h_CO_2^pure(T')+5.75
    h_H_2O^pure(T')}k_BT'^2dT'
    - [k_BTln{x^eq_H_2O(T)}-k_BT_3ln{x^eq_H_2O(T_3)}]




Note that Eq. (<ref>) is consistent with the view of Kashchiev and Firoozabadi <cit.> and with the reaction given by Eq. (<ref>) in which the hydrate is a new compound formed from x_occ molecules of CO_2 and 5.75
 molecules of water when the hydrate with occupancy x_occ.

Fig. <ref> shows the comparison between the driving force for nucleation obtained using the route 1 and a cutoff distance of r_c=1.0 nm when the hydrate is fully occupied and when only half of the small or D cages are occupied by CO_2 molecules. It is important to remark here that we are calculating driving forces for nucleation per cage of hydrate. This allows to compare Δμ_N^EC for both hydrates at the same conditions since the number of water molecules that form both solids is the same. Particularly, it is possible to know if a hydrate fully occupied is thermodynamically more stable than a hydrate partially occupied (85.5%) when both are formed from an aqueous solution of CO_2 at fixed conditions of pressure and temperature.

As can be seen, Δμ_N^EC is similar in both cases for low supercoolings. However, as the supercooling increases, the differences between both values increases. Particularly, Δμ_N^EC becomes less negative (driving force for nucleation is lower) when the hydrate is partially occupied than when the hydrate is fully occupied. This means that the fully occupied hydrate is more stable, from the thermodynamic point of view, than the hydrate with occupancy of 87.5% since the driving force for nucleation is higher. However it remains to be studied in the future if an occupancy between 0.875 and 1 could be more stable than the fully occupied hydrate.





Finally, it is also interesting to compare the driving force for the CO_2 hydrate with that corresponding to ice Ih obtained previously by Espinosa and coworkers. <cit.> They have obtained a value of Δμ^EC_N≈ -0.29 k_BT per water molecule considering a supercooling of 35 K at 1 bar with r_c=0.9 nm. For the case of the CO_2 hydrate, the driving force is -2.1 k_BT per unit of hydrate molecule at the same supercooling and 400 bar. This can be transformed to give -0.36 k_BT per unit of water molecule. Comparing both driving forces for nucleation, the nucleation of the CO_2 hydrate is more favorable than that of the ice Ih.

The previous result is true if the other factors that affect the nucleation rate are equal or similar, particularly the interfacial free energy. The ice Ih–water interfacial free energy for the TIP4/Ice model has been recently determined by Espinosa and collaborators <cit.> using the Mold Integration methodology. They found a value for the interfacial energy of 29.8(8) mJ/m. This value corresponds to the interfacial free energy averaged over all crystal orientations. Algaba et al. <cit.> and Zerón et al. <cit.> have also determined the CO_2 hydrate–water interfacial free energy using the Mold Integration-Host and Mold Integration-Guest. Particularly, they have found 29(2) mJ/m^2 and 30(2) mJ/m^2, respectively at 400 bar and 287 K. According to this, the nucleation of the CO_2 hydrate is more favourable than that of ice Ih since the interfacial energies of the ice and the CO_2 hydrate are similar.






§ CONCLUSIONS


In this work, we have studied the solubility of CO_2 in aqueous solutions when they are in contact via planar interfaces with a CO_2-rich liquid phase and with the hydrate phase at 400 bar using molecular dynamics computer simulations. We have also estimated the driving force for the nucleation of the CO_2 hydrate using four different routes. These properties are key to understanding, from a thermodynamics point of view, the parameters that control the nucleation of CO_2 hydrates. Water is described using the TIP4P/Ice water model and CO_2 using the TraPPE model. The unlike dispersive interactions between water and CO_2 are taken into account using the approach proposed by us several years ago. This selection allows to describe very accurately, not only the dissociation temperature of the hydrate at the pressure considered in this work, but also the CO_2 hydrate–water interfacial free energy. Calculations of solubilities have been carried out using the direct coexistence technique between two phases. Additional simulations of the pure systems, at several temperatures, have also been performed to calculate the driving force for nucleation along the isobar considered in this work.

We have analyzed the aqueous solution of CO_2 when it is in contact with the liquid phase (pure CO_2) and with the hydrate using two different values of the cutoff associated with the dispersive interactions. From this information, we have obtained the solubility of CO_2 in water when the solution is in contact with the CO_2 liquid phase. The solubility of CO_2 decreases with temperature, in a similar way to that of methane. However, the solubility of CO_2 is one order of magnitude larger than that of methane. We also observed an important effect of the long-range dispersive interactions in the solubility curve along the isobar of 400 bar. The solubility of methane in water is also affected by these contributions but their effect is smaller. It is interesting to remark that corrections due to the long-range dispersive interactions affect in a different way both systems. Whereas the solubility of CO_2 increases with the cutoff distance, in the case of methane it decreases. This is probably an effect due to the CO_2-CO_2 and CO_2-water electrostatic interactions. We have also studied the solubility of CO_2 in the aqueous solution when it is in contact with the hydrate and analyzed its interfacial structure. This magnitude increases with the temperature, as it happens with the solubility of methane in water. Contrary to what happens when the aqueous solution is in contact with the CO_2 liquid phase, the variation of the cutoff distance due to the long-range dispersive interaction has no effect on the solubility. This behavior has been also observed in our previous study dealing with the solubility of methane in water.

The dissociation temperature of the CO_2 hydrate (T_3), at 400 bar, can be evaluated from the intersection of the two solubility curves obtained in this work. This intersection is possible because the formation of the hydrate phase, at T<T_3, and the formation of the CO_2 liquid phase, at T>T_3, are activated processes. This means that there exists metastability below and above the dissociation temperature of the hydrate at 400 bar, and because of this one can find the intersection between the two solubility curves. The temperature at which this occurs is the T_3 of the hydrate at the fixed pressure. From this analysis we find that the dissociation temperature of the hydrate is located at 290(2) K when the cutoff distance for dispersive interactions is equal to 1.0 nm. This is in good agreement (within the error bars) with our previous estimation of T_3 obtained from direct coexistence simulations and using the same cutoff distance, 287(2) K. If the cutoff distance is larger (1.9 nm), the T_3 is located at 292(2) K. Although the value obtained in this work for a cutoff distance of 1.0 nm compares well with our previous estimation, 287(2) K (within the error bars), it is possible that finite-size effects produce a shift of T_3 towards higher temperatures, as well as the value used to account for the long-range dispersive interactions.

We also estimate the driving force for nucleation of the CO_2 hydrate. Particularly, we have calculated Δμ_N using the three routes proposed in our previous paper (routes 1-3). <cit.> Since the solubility of CO_2 in water is higher that that of methane by one order of magnitude, we have proposed a novel and alternative route based on the use of the solubility curve of CO_2 with the hydrate. This new route (which we refer to as route 4 in this paper) considers rigorously the non-ideality of the aqueous solution of CO_2 and provides reliable results for Δμ_N. Routes 1, 3, and 4 provide similar values of the driving force for nucleation of the CO_2 hydrate in a wide range of supercoolings. Unfortunately, the route 2 can not be used for CO_2 hydrates due to the non-ideality of the water + CO_2 mixture at the conditions considered.

Finally, we have also analyzed the effect of the cutoff distance due to dispersive interactions and the occupancy of the cages on the driving force for nucleation of the CO_2 hydrate. In both cases, there is a non-negligible effect on the driving force for nucleation. Particularly, the driving force for nucleation increases when the cutoff distances increases and when the occupancy of the small or D cages of the hydrate increases from 87.5% of occupancy to 100%.


The driving force for nucleation of the CO_2 hydrate obtained in this work using the novel route 4 is Δμ_N≈ -0.36 k_BT (per water molecule) at 400 bar and a supercooling of 35 K. This value lies below (is more negative) the driving force for nucleation of ice Ih at 1 bar and the same supercooling, Δμ_N≈ -0.29 k_BT (also per water molecule). Since other factors that affect to the nucleation rate are similar (i.e., the interfacial free energy), our results indicates that CO_2 hydrates nucleate more easily than ice Ih. However, should be taken into account that this comparison is not straightforward since CO_2 is needed for the formation of the CO_2 hydrate, but not for the ice.





§ CONFLICTS OF INTEREST

The authors have no conflicts to disclose.



§ ACKNOWLEDGEMENTS

This work was finnanced by Ministerio de Ciencia e Innovación (Grant No. PID2021-125081NB-I00), Junta de Andalucía (P20-00363), and Universidad de Huelva (P.O. FEDER UHU-1255522 and FEDER-UHU-202034), all four cofinanced by EU FEDER funds. We also acknowledge the Centro de Supercomputación de Galicia (CESGA, Santiago de Compostela, Spain) for providing access to computing facilities. The authors also acknowledge Project No. PID2019-105898GB-C21 of the Ministerio de Educación y Cultura. We also acknowledge access to supercomputer time from RES from project FI-2022-1-0019. J. G. acknowledges the nacional support from Gdansk University of Technology by the DEC-09/2021/IDUB/II.1/AMERICIUM/ZD grant under the AMERICIUM - “Excellence Initiative - Research University” program. Part of the computations were carried out at the Centre of Informatics Tricity Academic Supercomputer & Network. The research was supported in part by PL-Grid Infrastructure.



§ DATA AVAILABILITY


The data that support the findings of this study are available within the article.





