AN EXAMINATION OF THE PREDICTION OF HYDRATE FORMATION CONDITIONS IN THE PRESENCE OF THERMODYNAMIC INHIBITORS

Gas hydrates are crystalline compounds, solid structures where water traps small guest molecules, typically light gases, in cages formed by hydrogen bonds. They are notorious for causing problems in oil and gas production, transportation and processing. Gas hydrates may form at pressures and temperatures commonly found in natural gas and oil production pipelines, thus causing partial or complete pipe blockages. In order to inhibit hydrate formation, chemicals such as alcohols (e.g., ethanol, methanol, mono-ethylene glycol) and salts (sodium, magnesium or potassium chloride) are injected into the produced stream. The purpose of this work is to briefly review the literature on hydrate formation in mixtures containing light gases (hydrocarbons and carbon dioxide) and water in the presence of thermodynamic inhibitors. Four calculation methods to predict hydrate formation in those systems were examined and compared. Three commercial packages (Multiflash®, PVTSim® and CSMGem) and a hydrate prediction routine in Fortran90 using the van der Waals and Platteeuw theory and the Peng-Robinson equation of state were tested. Predictions given by the four methods were compared to independent experimental data from the literature. In general, the four methods were found to be reasonably accurate. CSMGem and Multiflash® showed the best results.


INTRODUCTION
Gas hydrates are ice-like, crystalline molecular structures which consist of two groups of components: host and guest molecules.The host molecule is water, whose hydrogen bonds form structures that trap guest molecules under suitable pressure and temperature conditions (Sloan and Koh, 2008).Hydrates are one of the main flow assurance problems faced by the oil and gas industry.Once hydrates form, they may agglomerate and produce plugs that either impair or interrupt the flow in the pipelines, bringing production to a halt in the latter case.This can occur in gas, gas/condensate, and oil production lines.
Amongst the techniques aimed at preventing hydrate formation, removing the free and dissolved water from the fluids to be transported is the most reliable one, albeit unfeasible in many situations.Keeping high temperatures and low pressures so that hydrates do not form is another strategy.However, in most offshore operations, hydrate formation is controlled by the injection of a thermodynamic hydrate inhibitor, such as alcohols (e.g., ethanol, methanol, mono-ethylene glycol, etc) and salts (e.g., sodium, potassium or calcium chloride).The prevention of hydrate formation and the safe removal of hydrate plugs incur high costs, raise safety concerns and are well-known deepwater flow assurance issues.
Amongst the common gases produced in the petroleum industry, methane, ethane, isobutane, propane, nitrogen, hydrogen sulfide and carbon dioxide are all hydrate formers (Carroll, 2004).It is well-known that carbon dioxide (CO 2 ) is one of the most common non-hydrocarbon gases found in petroleum reservoirs.In addition, most of the Pre-Salt fields in Brazil hold a significant amount of CO 2 , as high as 80% (Melo et al., 2011;Jahn et al., 2012).The pressure and temperature conditions in deepwater petroleum production (as that found in Pre-Salt fields) are very favourable to hydrate formation.Thus, both experimental and theoretical studies on hydrate formation, especially when high carbon dioxide contents are present, are important to the oil and gas industry.
Knowing and understanding the role of the pressures and temperatures in gas hydrate formation is the first step in controlling hydrate formation.A number of commercial simulators can predict the pressure/temperature conditions that favour hydrate formation.These simulators may use different thermodynamic models, especially when it comes to the calculation of the chemical potential of the aqueous and vapour phases, the effect of hydrocarbons of higher molecular weights, the impact of hydrate inhibitors, and the equation of state to obtain the fugacities of each component in each phase.A solid understanding of the models used in the available simulators is required to assess the quality of the results given by those simulators as far as the hydrate formation is concerned.
Only a few studies in the literature present a comparative assessment of models/simulators for hydrate prediction.Carroll (2004) examined and compared seven methods for predicting hydrate formation in mixtures containing hydrogen sulfide.For the set of data examined (125 experimental points), the computer methods CSMHYD and Equi-Phase showed to be reasonably accurate.The authors concluded that these two methods were able to predict the hydrate formation temperature to within 1.7°C in 90% of the cases.Ballard and Sloan (2004) compared the predictions from CSMGem, the hydrate program from the Colorado School of Mines, with CSMHYD (predecessor to CSMGem), and three commercial hydrate prediction programs (DBRHydrate, Multiflash ® , and PVTsim ® ).The comparisons were split into two categories for hydrate formation temperatures and pressures, uninhibited and thermodynamically inhibited systems.The authors concluded that all simulators fared favourably, giving results within about 1°C.
In this work, motivated by the high concentration of carbon dioxide in the Brazilian Pre-Salt fields, we compared the performance of four simulators to predict the hydrate phase equilibrium in systems containing carbon dioxide and light hydrocarbons, including systems with alcohols and salts (thermodynamic hydrate inhibitors).Due to the importance of high carbon dioxide content systems and the need to have reliable predictions of hydrate phase equilibria with and without thermodynamic inhibitors, this study gives a critical evaluation of the simulators.
The objective of this work is to briefly review the literature on hydrate formation in mixtures containing carbon dioxide and light hydrocarbons in the presence of thermodynamic inhibitors.Four simulators that predict hydrate formation for these systems are examined and compared here: three commercial products (Multiflash ® , PVTsim ® and CSMGem) and an in-house hydrate prediction program that uses the van der Waals and Platteeuw model (van der Waals and Platteeuw, 1959) coupled with the Peng-Robinson equation of state (EoS) (Peng and Robinson, 1976).The purpose is to compare hydrate prediction programs for the hydrate formation pressure at inhibited and uninhibited conditions, focusing mainly on systems containing carbon dioxide.In this work, the evaluated experimental data involved pure CO 2 hydrates, CO 2 -rich (>90 mol %) gas mixtures and also gas mixtures where CO 2 is considered a contaminant (<3 mol %), such as in natural gas compositions.The hydrate thermodynamic inhibitors investigated were alcohols and salts.

CALCULATION METHODS
The hydrate prediction program comparisons are divided into two categories for hydrate formation pressure: (1) uninhibited systems, and (2) thermodynamically inhibited systems.The first three simulators (CSMGem, Multiflash ® and PVTsim ® ) are based on rigorous thermodynamic models.They simulate hydrate formation conditions for gas and oil mixtures and can deal with the most commonly used hydrate thermodynamic inhibitors (methanol, ethanol, glycols, and salts).If the aqueous phase contains salts, the influence of the salts on the hydrate formation conditions will also be taken into account.The programs calculate the equilibrium conditions, amounts and types of hydrate structures (I, II and H), and composition of the phases present (i.e., hydrocarbon, aqueous, and solids phases).Details of the hydrate model from NUEM/UTFPR are presented in the next subsection.
NUEM/UTFPR Hydrate Model Kakitani (2014) implemented a model to predict the hydrate equilibrium formation.The basic model assumption relates to the chemical potential of water.At equilibrium, the chemical potential must be the same in all existing phases, that is, (1) where superscript H refer to the hydrate phase, L to the liquid phase, V to the vapour phase, and the subscript w to the water.
The chemical potential of the water in the liquid phase is obtained from classical thermodynamics as where f is the fugacity and the superscript 0 refers to a reference state for water.
The water chemical potential in the hydrate phase was obtained from statistical thermodynamics, namely from the van der Waals and Platteeuw (1959) model, as: (3) where the superscript β refers to a hypothetical phase in a crystalline configuration with empty water cavities, v i is the number of cavities of type i per water molecule in the hydrate structure and Y ki is the fractional occupancy of component k in cavity type i.In this work, the water chemical potential in the vapour phase was not considered since it was assumed that the solubility of water in the vapour is negligible.
Substituting Eqs. ( 2) and ( 3) in ( 1), and applying the Gibbs-Duhem equation, the equation in terms of measurable variables (p, T) becomes: where Δμ 0 , ΔH 0 and ΔV 0 are the differences in chemical potentials, molar enthalpies and molar volumes, respectively, between water in empty cavities and the reference state.These properties are tabulated in Table 1.The heat capacity difference between the empty hydrate lattice and pure liquid water phase is given by Δc p = -37.32+0.179(T-T0 ) according to Chapoy et al. (2012).
Equation ( 4) must be solved so as to obtain the state conditions for the hydrate equilibrium formation for a specific gas (guest).As Eq. ( 4) is implicit in pressure (p), an iterative process to calculate the equilibrium pressure at any temperature (T) becomes necessary.This equation was implemented as a Fortran90 code and solved by means of the secant method.
Fugacities of the gaseous components are calculated by the Peng and Robinson (1976) equation of state.Fugacities of water in the liquid phase are calculated by ( 5) From the definition of activity coefficient in a liquid mixture, (6) where x w is the mole fraction of the free water available to form hydrates.
Thus, the water fugacity in a liquid mixture in the presence of an inhibitor is: If there is no addition of inhibitors, that is, for an ideal solution, the activity coefficient is unity and the water fugacity is proportional to the mole fraction of the free water.If the solution is non-ideal, as is the case for mixtures with alcohol and salts, the activity coefficient must be calculated.Moreover, the gas solubility in water (x g = 1 − x w ) is usually considered to have no influence on the activity coefficient of water.As the gas solubility in the liquid water increases, the mole fraction of the water decreases, thus decreasing its activity.
In this work, two types of thermodynamic inhibitors were considered: alcohols and salts.Salts dissolve in water and form an electrolyte solution, which is different from the alcohol and water solution.This fact impacts the water activity coefficient calculation and in this work two activity coefficient models were used depending on the type of inhibitor.For alcohols, the water activity coefficient was calculated using the UNIQUAC model (Abrams and Prausnitz, 1975).For electrolytes, the water x f γ = activity coefficient was calculated using the Debye-Hückel model (Sander et al., 1986).Sloan and Koh (2008) proposed the use of the Krichevsky and Kasaranovsky (1935) correlation to calculate the mole fraction of the water.This correlation is used in this work.It follows, (8) where V is the gas molar volume in infinite dilution and H is the Henry constant that can be calculated by the following equation: In this equation, 0,1,2,3...

H
are the specific constants of each gas or guest.These particular constants were obtained from Sloan and Koh (2008).
The fractional occupancy of component k in the cavity type i, (Y ki ), is described by: (10) where C ki is the Langmuir constant for the gas component k in cavity type i and f k is the fugacity of the component k in the vapour phase.Parrish and Prausnitz (1972) improved the van der Waals and Platteeuw (1959) and McKoy and Sinanoglu (1963) models.They used the following expression to calculate the Langmuir constants: (11) where ω(r) is the intermolecular potential between the guest and the host, k is the Boltzmann constant, and r is the radial position.Parrish and Prausnitz (1972) and McKoy and Sinanoglu (1963) used the Kihara potential to account for the intermolecular potential.Munck et al. (1988) and Parrish and Prausnitz (1972) proposed the following expression to calculate the Langmuir constants: (12) where A ki and B ki are tabulated values adjusted to experimental data.
In this study, the Langmuir constants were calculated in two ways: using the expression proposed by Munck et al. (1988) and Parrish and Prausnitz (1972) (Eq.12) and using the Kihara potential (Eq.11) (McKoy and Sinanoglu, 1963).Table 2 shows the Kihara parameters (a, σ, ɛ/k) and parameters A and B (Munck correlation) for the calculation of Langmuir constants used in this work.
For further details about the model and solution algorithm in the implementation, see Kakitani (2014).

DATA ANALYSIS
Predictions given by the four methods were compared against experimental data from the literature.Table 3 shows the investigated systems and the related literature.All the experimental points investigated are below the upper quadruple point, that is, they are in the aqueous liquid-hydrate-vapour (L w -H-V) equilibrium.No data above the upper quadruple point were considered.In this study, a total of 225 experimental points are examined.
The four methods discussed above were used to calculate the hydrate equilibrium pressure, at a fixed temperature, point-by-point and the results are presented in terms of a percentage average absolute deviation, AARD, (13) where n is the number of points, p i exp is the experimental pressure and p i cal is the calculated pressure for point i.

RESULTS AND DISCUSSION
Each of the calculation methods discussed above was used to estimate the hydrate equilibrium pressure for each of the data sets.Figure 1 shows the AARDs of each calculation method for uninhibited systems.Figure 1(a) is for a single guest and Figure 1(b) is for more than one guest.Figure 1(a) and 1(b) suggest that one can expect the incipient hydrate pressure formation to be predicted within 11% and 20% accuracy for a single guest and for more than one guest, respectively.However, from Fig. 1(b), it is apparent that more reliable thermodynamic data are needed for systems with more than one guest.The AARDs for CSMGem, Multiflash ® and PVTsim ® are about 2%, 3% and 3.5% for a single guest, respectively, and 13.5%, 11% and 10% for systems with more than one guest, respectively.It can be observed that, for systems with a single guest and with more than one guest, the NUEM-UTFPR model that uses the Munck correlation performed better than the model that uses the Kihara potential.(Sloan and Koh, 2008).(Sloan and Koh, 2008) and Munck (Munck et al. 1988) parameters for the calculation of Langmuir constants.Figure 2 is a parity plot for the various prediction methods using uninhibited hydrate data.This graph is a plot of the predicted pressure as a function of experimental hydrate equilibrium pressure.If the prediction method were a perfect fit of the experimental data then all of the points would lie on the x=y line.Also plotted on the graph are error bands that deviate from the experimental pressure by ±10%.In general, this parity plot does not reveal much more than the error graphs presented earlier.However, the plot demonstrated that these methods generally tend to under-predict the hydrate equilibrium pressure, with the exception of the Fortran routine using the Kihara potential.In this plot the majority of the points are below the x = y line.It also can be observed that even at higher pressures the inaccuracies in the three commercial packages tend to be inside the error range of ±10%.

Gas
Figure 3 compares the prediction for a system containing thermodynamic inhibitors for a single-guest hydrate.Figure 3(a) and 3(b) are for inhibited systems with salts and alcohols, respectively.As expected, the AARDs for the inhibited systems are larger than those for uninhibited systems for all calculation methods.From these figures, one can expect commercial packages to predict the salt-and alcohol-inhibited incipient pressure within 9% approximately.Comparisons for the NUEM-UTFPR model for the incipient pressure are within 19% approximately.For the inhibited systems the model using Munck correlation performed slightly better than the Kihara potential.The type and amounts of salts and alcohols were included in Table 3.
Figure 4 is the incipient pressure for inhibited hydrate data for systems with more than one guest.One might expect commercial packages to predict inhibited incipient pressure within 10% approximately.Multiflash ® showed the best results.Figure 4 also shows that the inaccuracies for gas mixtures are higher than for a single guest in inhibited systems.For the NUEM-UTFPR model, the incipient pressure is predicted within about 20%.Once again, the model using Munck correlation performed slightly better than the Kihara potential.
Figure 5 shows the parity plot for the predictions for the inhibited hydrate equilibrium data.The plot shows that the models considered tend to over-predict the hydrate equilibrium pressure, with the exception of the NUEM-UTFPR model using the Munck correlation.The inaccuracies in the three commercial packages tend to be inside the error ranges of ±10% even at higher pressures.The NUEM-UTFPR model using the Kihara potential fails for pressures above 20,000 kPa.
In a general, CSMGem showed the best results for uninhibited one-guest systems, followed by Multiflash ® and PVTsim ® .However, for uninhibited more than one guest systems, PVTsim ® showed the best results followed by Multiflash ® and CSMGem.For inhibited systems the AARD of the three commercial packages are very similar when analysing one-guest systems and Multiflash ® worked better with systems with more than one guest.The NUEM-UTFPR model was less accurate than the three commercial software evaluated.
It must be pointed out that all calculation methods used the van der Waals and Platteeuw (1959) model for the hydrate phase.However, the fluid phases are treated by different equations of state.CSMGem uses the Soave-Redlich-Kwong (SRK) Eos, Multiflash ® uses the Cubic- Plus-Association (CPA) EoS, and PVTsim ® uses the Peng-Robinson (PR) EoS.This information, together with the type of activity coefficient equation used in the models could assist in understanding the differences between the calculation methods.Unfortunately, the commercial packages are not explicit with regard to the activity coefficient models considered.Certain activity coefficient models are more suitable to deal with specific   inhibitors like salts or alcohols than others.As for the equations of state, it is well known that cubic EoS are more conservative and general, despite the good results for hydrocarbon fluid phase equilibria.The CPA EoS combines the classical simple SRK equation with an association term (Kontogeorgis et al., 1996).Thus, this EoS has a "physical" part from the cubic EoS and a "chemical" part from the Association Theory.CPA has attracted the interest of the oil and gas industry, not just because of successful simulations over extended temperature and pressure ranges, but also because of the capability of the model to predict multicomponent, multiphase equilibria.In hydrate models, CPA represents successfully the partitioning of the inhibitor between vapour and aqueous phases.
The NUEM-UTFPR model to predict hydrate equilibrium formation presented the largest deviations.However, for engineering purpose, it could be a starting tool for estimating hydrate equilibrium predictions.

CONCLUSIONS
After reviewing some available carbon dioxide hydrate equilibrium data in the literature, four calculation methods to predict hydrate formation were examined and compared.The results presented here provide insight about the methods chosen.For the set of data examined in this study, observing the ranges of pressure, temperature and compositions investigated, totalling 225 experimental points, the simulators CSMGem and Multiflash ® exhibited the best results overall.
The investigated methods for predicting hydrate conditions are reasonably accurate for the evaluated systems.The total absolute average deviations (AARD) in the predicted hydrate pressure are: 7.3%, 16.4%, 11.6%, 7.3% and 7.7% for CSMGem, NUEM-UTFPR model with Kihara potential, NUEM-UTFPR model with Munck correlation, Multiflash ® and PVTSim ® , respectively.Therefore, CSMGem shows the best results for uninhibited one-guest systems, followed by Multiflash ® and PVTsim ® ; for uninhibited more than one guest systems, PVTsim ® showed the best results.For inhibited systems the deviations of the three commercial packages are very similar in oneguest systems, and Multiflash ® worked better with gas mixtures.
The NUEM-UTFPR model presented the largest deviations amongst the four calculation methods evaluated.The results using the Munck correlation performed better than those obtained by using the Kihara potential in five of the six scenarios evaluated: uninhibited systems (single guest), uninhibited systems (more than one guest), inhibited systems (single guest for alcohols), inhibited systems (more than one guest), and in the last scenario with all hydrate equilibrium data evaluated.One limitation of this model could be the assumption that the solubility of water in the vapour is negligible.
All calculation methods performed better for the single guest hydrates than multi-guest hydrate systems.This fact may be related to the equation of state (EoS) used by each method and the challenges of dealing with multicomponent multiphase equilibria, including the binary interaction parameters used in the EoS in the case of mixtures.It can be added that the Kihara parameters for the calculation of Langmuir constants are usually optimized to pure (singlecomponent) hydrate experimental data.This optimization is essential to the success of the methods.
The calculation methods exhibited the best results with the uninhibited systems.The presence of inhibitors in these systems is represented by the activity coefficient models, which may have great influence on the results.The commercial packages are not explicit in the activity coefficient models used and, depending on the activity coefficient model and how its parameters were regressed, some may be more adequate to deal with specific inhibitors like salts or alcohols.
Based on this study, one can also conclude that is not possible to accomplish a complete comparison between the calculation methods with the carbon dioxide hydrate equilibrium data evaluated.Notwithstanding, the specific comparison realized in this work is valid in the particular ranges of pressure, temperature and inhibitor compositions informed.This work shows the limitations of the predictive tools considered, which is an important exercise to determine the data gaps and applicability of the predictions.A comparison of prediction tools used for this purpose and an in-house hydrate prediction program in a specific range of pressure, temperature and inhibitor compositions may contribute to control hydrate formation in these kinds of systems, since this kind of information suggests directions in the selection of appropriate methods for carbon dioxide containing systems.

Table 1 .
Physical properties of sI hydrate

Table 3 .
Experimental data used to compare the models.