Accuracy of different thermodynamic software packages in predicting hydrate dissociation conditions

Gas hydrate blockage is a longstanding problem in the oil and gas industry, responsible for the interruption of hydrocarbon ﬂow, posing economic and safety risks. Computational modeling methods are required to estimate the dissociation pressure of the original system and its shifting in the presence of thermodynamic inhibitors so as to avoid excessive inhibitors dosages, high operating costs and serious environmental impacts. On the other hand


Introduction
Clathrate hydrates, simply known as gas hydrates, are solid solutions of hydrogen-bonded cavities of water molecules, in which small gas molecules are trapped.Typically, gas molecules smaller than 0.9 nm are encapsulated in the hydrate cavities, whose structural stability depends on the Van-der-Waals and London forces generated by the interactions between guest and host molecules [1] .Depending on the size of the guest molecules, three hydrate structures are known to form: structure I (sI), II (sII) and H (sH) [2 , 3] .The required conditions for their formation are sufficient supply of gas and water under suitable pressure and temperature conditions [4 , 5] .
Gas hydrates have been studied for their potential as an energy source, as the amount of carbon trapped as natural gas in gas hydrate reserves has been estimated to be as high as twice the amount of carbon contained in all fossil fuels on Earth [6 , 7] .Due to their significantly high volumetric gas storage capacity [3] clathrate hydrates have been considered as a possible solution for several technological applications [7] .The use of hydrates as a molecular storage medium for hydrogen and natural gas has been extensively investigated [8][9][10][11][12] , and efforts to develop clathrate hydrate materials for hydrogen storage have proven that hydrogen can be stored either at high pressure conditions [10] or at lower pressures with the addition of promoters [11] .Since methane hydrates have an energy density equivalent to that of highly compressed https://doi.org/10.1016/j.ctta.2022 Coal Bed Methane natural gas, clathrate hydrates can also be an attractive solution for the transportation sector [12] .Hydrate-based gas separation technology aims to isolate CO 2 from flue gas streams [13] , thus revealing the possible application of hydrates to capture greenhouse gasses.The main advantages of this method compared to conventional adsorption and absorption CO 2 capture techniques is the lower energy demand and that it is more environmentally friendly [14 , 15] .Additional applications of hydrates include water treatment and desalination [16 , 17] refrigeration [18 , 19] and other separation processes [20 , 21] .
Although hydrates are a potential commodity and a promising technology for the future, they also pose a threat of severe plugging problems in the oil and gas industry [22] .With the progress of deepwater exploration and production of hydrocarbons, hydrate hazards have evolved, leading to inevitable challenges in drilling, completion, and operation of production pipelines [23] .Extreme pressure and temperature conditions prevailing in deepwater environments have increased the risk of hydrate formation, leading to additional well-control challenges [24] .When hydrate blockage forms in a pipeline, high pressure develops in the upstream section thus imposing pipeline fracture risks.Moreover, the pressure difference between the upstream and the downstream sides may force the plug to experience a projectile behavior.Both events can result in loss of production, damage to the pipeline, and safety hazards for the personnel and the environment [3] .
The standard industry practice is to mitigate hydrate formation risk by continuous injection of thermodynamic hydrate inhibitors (THIs), such as mono-ethylene glycol [25] .Although this constitutes an effective solution for shifting the hydrate equilibrium curve towards higher pressures [26] , this technique involves significant operational cost overheads to the pipeline network operators.As a result, any practical solution to reduce the amount of inhibitors, hence the operating costs, while keeping the facility risk-free is of absolute necessity.The simplest way to reduce the inhibitor quantity is to avoid overdosing by perfectly matching the inhibitor concentration with the required level of inhibition, plus a safety factor.For this task, accurate predictions of the hydrate formation pressure and temperature conditions must be used for the design of the inhibition strategy.
Fast and reliable prediction tools are utilized to determine the hydrate phase boundary of the gas/inhibitor mixture of interest.The efforts in this field can be divided into three main categories: 1) the K-value method 2) gas gravity diagrams 3) rigorous thermodynamic models.As far as the K-value method is concerned, Katz [27] and his coworker Wilcox [28] generated K-value diagrams for pure gasses describing the equilibrium of hydrocarbon, water and hydrates.Gas gravity diagrams are simple and straightforward representations of hydrate stability at various temperature and pressure conditions, where the only correlating parameter is the average molar mass of the mixture, thus obscuring the exact role of each gas component.This method provides low accuracy predictions of the stability conditions for hydrocarbon mixtures and becomes even worse when non-hydrocarbon gas components are considered [29] .
Rigorous thermodynamic modeling of hydrates firstly appeared when Van der Waals and Platteeuw (VdW-P) derived the basic statistical thermodynamic equation for gas hydrates [30] , based on the Lennard-Jones potential.Later, McKoy and Sinano ğlu 1963 [31] provided improved predictions of hydrate dissociation conditions, by assuming a linear or spherical Kihara potential [31] .Parrish and Prausnitz [32] combined the VdW-P model with the Kihara potential and extended the capability of the model to describe phase equilibria of simple and mixed gas hydrates [33] .The aforementioned methods are combined with fluid and pure solid phase models for the estimation of hydrate formation conditions in multiphase flash algorithms [34 , 35] .Mass balance calculations coupled with thermodynamic modeling and Gibbs energy global minimization are standard methods to predict hydrate formation conditions, equilibrating fluid phase compositions and properties.
A number of commercial PVT simulators are available to predict hydrate formation conditions using the rigorous thermodynamic approach.Multiflash, HydraFLASH, CSMGem and CSMHyd are strong solvers for hydrate phase equilibrium problems, which, however, differ in the thermodynamic models used to predict hydrate dissociation conditions.Some use the original VdW-P theory for the hydrate phase modeling, which does not consider the distortion of the lattice structure due to the introduction of the gas molecule, whereas others do so by utilizing the modified VdW-P theory.Moreover, some PVT simulators use the general thermodynamic - approach, where an activity coefficient model and an EoS are used to calculate the fugacity of water in the liquid and the vapor phase respectively, while others use the  - approach.Some software packages utilize the Cubic Plus Association (CPA) EoS, firstly presented by Kontogeorgis et al. [36] , to model all fluid phases while others use the classic cubic EoS for the hydrocarbon fluids.For the effect of electrolytes in the aqueous phase some methods adopt the Debye Hückel theory [37] and others combine the Debye Hückel term together with additional virial and Born terms [38 , 39] in order to include the contribution of short-range ions interactions to the interactions between ions and the surrounding medium.The fugacity models for the pure solid phases also vary as some adopt the solid freeze out model, whereas others introduce additional corrections (such as the Poynting correction factor) to improve the accuracy of the fugacity for the solid phases.Finally, the values of the tuned parameters for the hydrate model (i.e., Kihara potentials) vary between simulators, depending on the set of the experimental data which has been used by each developer to tune their model.
In this work, four commercial simulators for hydrate prediction are reviewed: Multiflash by KBC, HydraFLASH by Hydrafact, CSMGem and CSMHyd by E. Sloan.The comparative evaluation of the simulators is carried out by considering five thermodynamic models of them.To the best of our knowledge all commercial software packages are retuned on a regular basis which implies that the utilized thermodynamic models have already been tuned, verified and compared among classic simple systems with one component pure hydrate former or inhibited systems in low salt concentrations.As a result, their performance is quite similar and surely satisfactory as far as the predicted association pressure and temperature conditions are concerned [40][41][42][43] .To avoid unnecessarily repeating such results, our work is focused on comparing the prediction performance of the models against newer results on more complex multicomponent gas hydrate former systems and on highly inhibited systems.This way, the models are pushed to the limits of their prediction capability thus allowing for a fair comparison.The comparison results are divided into two main galleries, namely pure hydrate former mixtures and salts inhibited systems.Depending on the type of industrial application, the pure systems gallery is further divided into two collections: flow assurance and technical application systems.All thermodynamic approaches are compared against independent experimental data which has been collected from the recent literature published between 2015 and 2019.
The rest of the paper is organized as follows.Section 2 focuses on the modeling approaches available to estimate hydrate dissociation pressure, the rigorous thermodynamics models and the solution algorithms.Section 3 describes the characteristics of each commercial software package, followed by thorough analytics of the data collected and utilized.Section 4 presents the results of each software against all galleries and collections, followed by a detailed analysis of their performance.The paper ends with the conclusions obtained in this work.

Hydrate formation modeling approaches
The rigorous approach to estimate dissociation pressure (or temperature) requires that at dissociation conditions the fugacities of all components are equal across all phases present and hydrates appear as an incipient phase.Especially for the hydrate phase, only water potential is considered by the hydrate model.Therefore, we require that: where superscripts , ,  ,  correspond to the hydrate, liquid, vapor and ice phases (if any).To satisfy those conditions, dissociation pressure (or temperature) is iterated until the fugacity of each component, which are computed by a multiphase flash algorithm, equilibrate across the non-hydrate phases and, at the same time, the fugacity of water in the hydrate phase equals to that in the non-hydrate phases.A simpler, yet computationally efficient approach does not employ complex equilibrium calculations but rather iterates pressure (or temperature) until the rigorously computed chemical potential of water in the hydrate phase equals to the value obtained by correcting experimental values from the reference (  0 ,  0 ) to the operating conditions.Both approaches are discussed in more detail in the following paragraphs.Regardless of the solution algorithm utilized, both require thermodynamic data of all phases obtained by suitable EoS models.The thermodynamics of the hydrocarbon phases are described by classic cubic EoS models.When polar compounds such as water, salts and electrolytes coexist, then more complex equations of state can be employed, such as the Cubic Plus Association model (CPA), or modifications/additional models should be set on top of the classic cubic EoS.The chemical potential of water in the hydrates is obtained using the original Van der Waals-Platteeuw theory and its modifications.

Hydrate phase modeling
Hydrates thermodynamics modeling was firstly attempted by Van der Waals and Platteeuw in 1959 (VdW&P) who computed the chemical potential difference between water in the aqueous and in the hydrate phase Δ  −   =    −    .For that purpose, they conceptualized a theoretical phase state in which water can form an empty cage as the actual hydrates so that Δ  −   is now split into two steps so that , where descriptor  denotes the theoretical empty water cage which is considered as a metastable phase.For the first term, the chemical potential difference between the hypothetical hydrate phase and a pure water phase can be computed using rigorous thermodynamics, based on the known chemical potential difference between pure water and the empty hydrate lattice at some reference conditions  0 ,  0 (such as the triple point of water).For the Δ −H  term, they applied the physics of adsorption of gas molecules in the cage by utilizing Langmuir's theory, thus introducing the Langmuir constant for gas component  in cavity type  ,  , , which is a function of temperature  that accounts for the gas-water interaction in cavities and it is computed as a function of the cell potential of component  in cavity type  ,  (  ) , , usually based on the Kihara potential with a spherical core [44] .Later, Ballard and Sloan modified Van der Waals-Platteeuw by adding an activity coefficient for the water in hydrate,   , to account for the distortion of the lattice structure by guest molecules [36] .For more information, the reader is referred to [39] .
In a different approach, Parish and Prausnitz (P&P) [32] developed a method that utilizes available experimental data on the chemical potential difference of water.They developed experimentally originated expressions to estimate the difference Δ  −    by integrating Δ  −    (  0 ,  0 ) of some reference hydrate at some reference state  0 ,  0 to any desired conditions  ,  .For the reference values utilized, xenon hydrate is adopted as reference for the sI hydrate structures and temperatures below 0 °C, whereas methane hydrate is utilized for higher temperatures.Similarly, CBrClF 2 (bromochlorodifluoromethane) hydrate is utilized for sII hydrate structure and temperatures below 0 °C as well as for natural gas at temperatures above that.

Fluid phase modeling
Conventional cubic EoS coupled with association is nowadays the most pronounced way to model non hydrate phases in hydrate simulators.On top of the classic non-polar interactions modeled by SRK, typically used to describe non-aqueous, fluid phase components, the CPA EoS ( [36] ) further considers the contributions to the compressibility factor (eventually to the Helmholtz free energy) due to self-and crossassociation hydrogen bonding.Modeling of an aqueous phase becomes complicated when inhibitors are utilized as the interactions of salt, water and alcohol need to be considered in the fugacity of water and of the solutes.The Shock and Helgeson equation (1988) [45] with a coupled Bromley activity model (1973) [46] is an implicit fugacity model providing the chemical potential of aqueous phase species which is then manipulated to obtain the fugacity by defining a suitable reference on which the chemical potential is based.
Alternative approaches model the aqueous phase by adding extra terms to the total residual Helmholtz free energy   , that is    for the cubic EoS (e.g., SRK),   for the chemical association and   for contributions of electrolyte solutions.Some of them only consider    using an advanced version of SRK/RKS combined to the NRTL mixing rule whereas the salt effect on hydrate formation is modeled using a pseudo equivalent-NaCl salt component approach.Other models account for all contribution terms to the Helmholtz energy and further expand   using various extra models.Those terms include   − ≋ to describe the electrostatic ion-ion interactions,    to describe the short-range interactions between ions and   to account for the interactions between ions and the surrounding medium.For more information on the development of Debye-Hückel, virial terms and Born's term the reader is directed to [37-39 , 47] .

Solid phase modeling
Two models are used in the tested algorithms to consider solid phases appearing when a liquid, such as water, freezes out, together with the hydrates.The first involves the use of the solid freeze-out model, which estimates the thermodynamic properties of fixed composition solid phases by freezing one or more components of the fluid mixture.The second model adopts the Poynting correction to model fugacity of ice which utilizes the fugacity coefficient of water in the vapor phase.

Estimation of hydrate dissociation conditions
At dissociation conditions: 1. equilibrium phase compositions honor the mass balance 2. each component exhibits equal chemical potential across all equilibrium phases 3. the total Gibbs energy is at its global minimum 4. hydrate phase appears as an incremental phase To attack the numerical problem, the feed is split at current dissociation pressure (or temperature) estimate using any powerful multiphase flash algorithm to get the fraction, composition and fugacity of each non-hydrate phase in equilibrium.Subsequently, the water fugacity in the hydrate phase is estimated at current conditions using the VdW-P theory or its modified version and it is compared to that of the coexisting phases.If all water fugacities are equal, the algorithm is terminated, otherwise the dissociation pressure (or temperature) is updated to eliminate the fugacities difference.The classic cost function of Michelsen [35] , exhibiting very nice convexity properties which facilitate convergence, and the one developed by Gupta et al. [34] are the two multiphase flash algorithms implemented in the software packages examined in this work.
Alternatively, the Parrish and Prausnitz (P&P) approach takes advantage of the fact that at the dissociation pressure the experimentally determined Δ  −    ( Section 2.1 ) must equal to the rigorously calculated value.By incorporating an activity coefficient model of the form Δ  −  =  ⋅  ln (    (  ,  ) ) , they were able to rigorously estimate

Evaluated thermodynamic approaches/models
The hydrate modeling approaches, and the fluid phase thermodynamic models discussed above are implemented in the four commercial software packages for phase equilibrium evaluated in this work, as shown in Table 1 .

HydraFLASH CPA (HF CPA)
For a given water/hydrocarbon/inhibitor mixture, HydraFLASH provides the most stable hydrate structure (sI, sII or sH) in equilibrium and describes the inhibition effect of salts/alcohols with an option to estimate the minimum inhibitor injection rate to prevent hydrates.It utilizes Michelsen's multiphase flash algorithm to predict hydrate formation conditions and the original VdW-P theory coupled with the CPA EoS to implement a fugacity-based model.For the non-associative term, it supports the SRK and the PR cubic EOS with various options for the alpha function in the attractive term including the formulation of Twu et al. [48 , 49] , the Mathias-Copeman one, as well as the classical Soave expression [50] .In addition, HydraFLASH supports the classic Van der Waals and the Huron-Vidal mixing rules [51] .
To ensure a fair comparison with the other CPA models tested SRK-EoS was selected as the non-associative part of the CPA which is further supported by the literature and by previous experience on satisfactory results in modeling hydrates systems [52] , especially when inhibition effect is considered [39 , 40] .The Soave model was selected for the attractive term, as it has been shown to lead to a consistent behavior in the prediction of both saturation properties as well as enthalpies and heat capacities [53] .The Huron-Vidal mixing rule was utilized due to its superiority in modeling the energetic effect and the asymmetry of simple and complex mixtures [44 , 45 , 54] .HydraFLASH models electrolyte solutions by adding only the Debye Huckel term to the total residual Helmholtz free energy [55] and utilizes the Poynting factor to correct the fugacity of water and to accurately model that of the ice phase.

Multiflash CPA (MF CPA)
Multiflash, by KBC, offers two thermodynamic models for hydrate systems, namely the CPA based one (labeled as MF CPA in this work) and the advanced version of RKS with an NRTL type mixing rule (labeled as MF RKSA in this work).Both models are capable of predicting hydrate dissociation/formation conditions as well as the amount of hydrates formed, the type of gas hydrate in equilibrium (sI, sII and sH) and the phase transitions among the fluid and hydrate phases.Multiflash supports modeling of inhibition effects as well as the estimation of the required inhibitor injection rate.
Similarly to HF CPA, it is a fugacity-based model utilizing Michelsen's multiphase flash algorithm combined to the original VdW-P theory and the CPA EoS.MF CPA embeds the SRK EoS as the nonassociative term option, with a polynomial function to describe the alpha parameter and standard Van der Waals mixing rules, combined to the generic numerical formulation of CPA.To model the electrochemical effects of aqueous salt solutions the Debye-Huckel method is used followed by virial corrections and a Born term [38 , 39 , 47] .MF CPA model uses the freeze-out model to calculate phase equilibrium and properties of other solid phases that may coexist with hydrate (e.g., ice).

Multiflash Redlich Kwong Soave advanced (MF RKSA)
To enhance the fugacity calculation of complex aqueous phases (typically containing electrolytes, alcohols), MF RKSA utilizes Michelsen's multiphase flash algorithm, this time coupled with the Advanced Redlich Kwong Soave (RKSA) EoS and the NRTL mixing rule type.MF RKSA considers the residual Helmholtz free energy of the system as a function of the SRK term with no addition of association or electrolyte ones as it is the case of MF CPA.The Advanced RKS considers an SRK attractive parameter which is determined by fitting vapor pressure to improve the saturation pressure prediction, it applies the Peneloux volume correction to improve the liquid density and utilizes the NRTL mixing rule apart from the Van der Waals classic one.Yet, the binary interaction coefficients of polar components   ≠   and ∝  = ∝  , hence when ∝  is set to zero, the NRTL mixing rule reduces to the conventional Van der Waals one with symmetric binary interaction parameters   =   .
For electrolytes the pseudo-salt component approach is employed assuming that salts are represented by a single pseudo NaCl-equivalent salt component, exhibiting its own physical properties and binary interaction parameters and it is assumed to behave as a non-volatile-heavy component in the aqueous phase, not to intervene with the hydrocarbon phase and to precipitate at the solubility limit.
For aqueous phase solutions containing various types of salts (e.g., MgCl 2 , CaCl 2 , NaCl etc.), user defined rules are assigned to the model to convert the mixed salt solution into a single equivalent pseudo-NaCl solution.Finally, similar to the MF CPA model, the solid freeze-out model is utilized to compute component fugacity of solid phases that coexist with hydrates (e.g.ice).

CSMGem (CSMGem)
This software package has been developed at the Center of Hydrate Research in Colorado School of Mines (CSM).CSMGem supports the modified VdW-P model in a fugacity-based model, coupled with SRK EoS to model hydrocarbon phases, the Shock -Helgeson equation and the Bromley activity model for the aqueous phase.It utilizes the flash algorithm of Gupta et al. [34] to perform multi-phase equilibria on the basis of Gibbs energy minimization.
CSMGem supports equilibrium calculations with up to 10 phases, depending on the type of components that are present in the feed, i.e. vapor and liquid hydrocarbon, aqueous, ice, hydrates (sI, sII and sH) as well as solid NaCl, KCl and CaCl 2 but not MgCl 2 .For a given feed, the multiphase algorithm initially accounts for all possible phases and subsequently their number is reduced along the iterations of the flash algorithm.CSMGem provides detailed information on the number of phases present, the hydrates structure, the cage occupancy at given pressure/temperature and the equilibrium phases composition.

CSMHyd (CSMHyd)
CSMHYD has been the predecessor of CSMGem predicting hydrate formation temperature/pressure of sI and sII structure, with or without methanol and/or salts, at three and four phase conditions.Eight salts are supported, namely NaCl, Na 2 So 4 , NaF, KBr, KCl, CaCl 2 , MgCl 2 and SrCl 2 .This software package utilizes the original VdW-P approach coupled with SRK EoS and a variation of the P&P method.The results output of the incipient hydrate phase contains all equilibrium phases (L w -H-V, I-H-V or L w -H-V-L hc ), hydrate structure type, hydrate dissociation/formation conditions, compositions of phases and hydrate cage fractional occupancy.Similarly to CSMGem, CSMHyd makes use of the SRK EoS coupled with a statistical thermodynamic model to predict hydrate dissociation conditions.Yet, this model is based on fugacity as it estimates hydrate pressure and temperature equilibrium conditions so as to match theoretical and experimental Δ values.

Experimental data
An independent database of hydrate dissociation points was generated by retrieving all recent experimental data that was published between 2015 and 2019.17 journal papers providing hydrate dissociation pressure and temperature for various stream compositions were utilized.The number of data retrieved from each paper varied from just a few (3 dissociation pressure data points) to dozens (58 data points in a single work), depending on the number of feed compositions that were handled in each work.Although regular operating conditions are of main practical concern, data on extreme pressures were also retrieved to allow for the evaluation of the commercial software packages against such extraordinary cases.
To organize the evaluation process, the collected data points were divided into two galleries, those of the uninhibited and of the inhibited systems.The data in the uninhibited systems gallery is further divided the flow assurance issues collection and to industrial applications one, ending up in four uninhibited systems.The inhibited gallery includes systems of methane hydrate formers with various salt types and a wide range of salt concentrations and is split into two collections based on whether one or more salts are considered in the aqueous phase thus ending up to five more systems.The organization of the retrieved data together with some basic analytics is shown in Table 2 .Detailed information on the test data points can be found in the Appendix.
System 1, CO 2 -CH 4 mixtures, refer to natural gas (NG) produced from oil and/or dry gas reservoirs containing CO 2 as an impurity component and their study has been associated with flow assurance problems in the gas industry and the transportation sector and with the exploitation of methane hydrate reservoirs in permafrost regions and subsea sediments under various operating conditions, focusing on the replacement of CH 4 with high-pressure CO 2 , which favors the long-term storage of carbon dioxide [56 , 57] .
The CO 2 -N 2 -CH 4 mixtures in System 2 are also associated with flow assurance issues.In addition, producing methane hydrate reservoirs by exposing then to CO 2 -N 2 mixtures has been investigated [58] and implemented in fields beneath the Alaskan permafrost [58] .In this case, N 2 acts as a carrier gas mitigating technical problems associated with the injection of dense CO 2 liquid into the formation, thus combining carbon storage and methane hydrate exploitation.System 3 components, i.e., CO 2 and N 2 , are the main two components of flue gas, hence related to carbon capture [59] .Moreover, such systems are related to the conversion of natural gas hydrate deposits to CO 2 for both future power generation and greenhouse gas control, by replacing CH 4 in the hydrates with a gas mixture of CO 2 and N 2 [58] .
CH 4 -N 2 mixtures in System 4 are related to flow assurance issues during production and transportation operations.Such mixtures are also found in raw coalbed methane (CBM) mining sites, where separation of CH 4 and N 2 by conventional methods is inapplicable as both components are in supercritical state.Several research groups have developed potential gas separation techniques that exploit the differences in equi- librium concentrations between the vapor and hydrate phases [3 , 60] , leading to hydrate-based pressure-controlled separation of CH 4 .Systems 5 to 9 correspond to inhibited mixtures found in subsea natural gas pipelines where chemical hydrate inhibitors such as methanol and mono-ethylene glycol are used in large quantities ( > 30 wt%) to anticipate the conditions which favor hydrate formation.The presence of salts in the produced formation water imposes the need to use robust and reliable simulation models to predict hydrate formation conditions in such an environment and over a wide range of operating temperatures and pressures.In this gallery, CH 4 mixtures with an aqueous phase containing a pure salt are considered.The salts studied are NaCl, KCl, MgCl 2 and CaCl 2 .

Results visualization
The simulation results were firstly evaluated by means of classic statistical measures.Subsequently, an alternative and more thorough investigation was run allowing the attribution of the obtained deviations to distinct factors.
In the first evaluation mode, the typical absolute and relative error definitions were utilized where    is the experimental pressure and   , is the calculated pressure for point  using thermodynamic model .It was decided to work on the pressure deviation itself rather than on its relative error as only the former is of use when designing pipelines and applying safety margins.In the second evaluation mode, the data proved to be quite sparse to allow for a decent interpretation by means of 3D plots.To tackle this problem, the experimental data points of each gallery/collection was split into clusters of equal temperature range, each 2 K wide, and the deviations in each cluster were averaged.The Mean Absolute Error (MAE) and the Standard Deviation of the Absolute Error (SDAE) temperature cluster  are calculated by where  (  ) is the number of data points in the cluster.Similarly, The Mean Absolute Relative Error (MARE) and the Standard Deviation of the Absolute Relative Error (SDARE) are given by All pressures are given in MPa's.

Results summary
In this section an overall evaluation of the performance of the tested software packages is attempted by utilizing the absolute error range, the average signed and absolute pressure error in each gallery/collection as well as their deviation.Those results are shown in the Tables which accompany the discussion on each evaluated system.Nevertheless, a set of Figures showing extra information such as the predictions error distributions, their best fit probability curve and important statistics (average pressure error , pressure standard deviation  and average absolute pressure error ||,  +  and  −  values) have also been included in the supplementary material section.

CO 2 -CH 4 mixtures
The dataset examined consists of 57 data points which correspond to dissociation pressures between 1.78 and 25 MPa ( Table 3 ).All software packages except CSMHyd achieved balanced predictions which were shown to shift by less than 0.2 MPa.The error distribution exhibited limited variance with standard deviation  less than 0.7 MPa and average absolute error || below 0.5 MPa.The errors obtained indicate that those methods can safely predict the dissociation pressure within an acceptable error range for flow assurance applications.
CMSHyd provided consistently lower predictions which pushed  to − 0.6 MPa and increased || to 0.7 MPa.Additionally, its worse prediction on the dataset was 5.76 MPa, compared to the 2.04-2.74MPa achieved by other packages.

CO 2 -N 2 -CH 4 mixtures
This dataset consists of 35 data points corresponding to dissociation pressures between 4.81 and 30.66 MPa.Descent performance was achieved by all software packages with HF CPA, MF CPA and MF RKSA providing slightly better results as the predictions were balanced with an  value of 0.5 MPa and  equal to 0.9 MPa ( Table 4 ).
CSMGem achieved slightly worse with = 0.7 and = 1.1 MPa respectively, thus leading to an increased || value of 1.1 MPa.Nevertheless, as the worse error | |  achieved is almost same for all methods ( Table 4 ), the dissociation pressure of such mixtures can be safely predicted by all software packages.

CO 2 -N 2 mixtures
77 data points were considered in this dataset corresponding to dissociation pressures between 2.05 and 55.11 MPa as shown in Table 5 .It is remarkable that all models based on the VdW-Platteeuw theory exhibit almost identical performance.Although the dissociation pressure is consistently underestimated (  = 1.7 MPa),  is of the order of 1.7 MPa rendering the models sufficiently reliable for most engineering applications.CSMHyd overestimates dissociation pressure by only 1 MPa but exhibits larger spread of the errors thus increasing  to 2.7 MPa.Despite the differences in their error distributions, all methods arrive to an || value of about 2.0 MPa.

CH 4 -N 2 mixtures
This dataset contained 23 data points with dissociation pressures varying in 5 to 25 MPa ( Table 6 ).Similarly to the CO 2 -N 2 mixtures, all models but CSMHyd exhibit remarkably close performance with an  value of 0.8 MPa indicating dissociation pressure overestimation and a limited  value of 0.9 MPa.MF CPA exhibits slightly improved performance compared to its competitors with  = 0.8 MPa.CSMHyd exhibits milder overestimation but, as it was the case in the previous system, the errors scatter is higher.Nevertheless, results show that CH 4 -N 2 systems can be handled with descent accuracy by all methods.

CH 4 -Saline aqueous phase mixtures (Pure salt: NaCl)
The dataset of NaCl to inhibited mixtures contains 56 points corresponding to pressures between 2.71 and 183.9 MPa ( Table 7 ) with the elevated pressures corresponding to systems of high NaCl concentration and high temperature.The CPA based methods exhibit similar performance with an offset of 5.3 and 3.7 MPa respectively, and slightly better performance for the MF CPA as  = 7.2 MPA compared to 10.0 MPA for the HF CPA.
CSMGem exhibits remarkably balanced and fine performance with  lowering down to only 4 MPa and |  | to 15.47 MPa, even for dissociation pressures of the order of 150 MPa or more.On the other hand, MF RKSA and CSMHyd exhibit very poor, highly biased performance with  = 16.4MPa or higher and || = 9.6 MPa or more, thus rendering them  as incapable of providing dissociation pressure predictions of descent accuracy.
4.1.6.CH 4 -Saline aqueous phase mixtures (Pure salt: KCl) 37 data points of KCl inhibited systems were tested in a pressure range of 4.02 to 180.06 MPa ( Table 8 ).The CPA based models exhibit almost identical performance with = 6.1 MPa and a very large deviation of 46.5 MPa.As it will be shown in the detailed analysis, this is due to their failure to perform in high salt concentration.MF RKSA exhibits similar || value (16 MPa) due to the combination of worse offset and improved error variance.CSMHyd exhibits errors spread over a very wide range, leading to a very high || value (23.8 MPa).CSMGem is remarkably accurate with a slight only offset, |  | of 24.31 MPa and  of only 5.1 MPa thus clearly rendering it as the most reliable tool to handle KCl inhibited systems.

CH 4 -Saline aqueous phase mixtures (Pure salt: MgCl 2 )
The performance of MF CPA, MF RKSA and CSMHyd against the 21 data points available is very poor with  being as high as 21.8, 35.8 and 30.4 MPa respectively ( Table 9 ).Additionally, CSMGem did not provide predictions as it does not support the MgCl 2 salt.On the other hand, HF CPA exhibited great accuracy with  of only 2.8 MPa and || of only 2.9 MPa, despite the very large values of the dissociation pressure.

CH 4 -Saline aqueous phase mixtures (Pure salt: CaCl 2 )
This dataset contained 29 data points which can be split into two groups ( Tables 10 ).The first group refers to CaCl 2 concentration below 23.8 wt% and regular temperatures lying above 273 K.The second group corresponds to 6 measurements run on a 32 wt% CaCl 2 sample at extremely low temperatures below 264 K.Most deviations lie close to the baseline except of the second group points which are spread towards much higher deviation values.As will be shown in the detailed results section, most evaluated methods performed fine in the first data points group while they exhibited significant deviations in the second one.
Despite this discrepancy, the CPA based methods perform slightly better compared to the CMSGem and CMSHyd ones with an || value of 10.0 MPa.MFRKSA cannot handle the second group measurements at all, the deviations of which drive its bias to extremely low values (32.5 MPa) and || to a very high value of 32.8 MPa.

CH 4 -Saline aqueous phase mixtures (including salt mixtures)
36 data points with salt mixtures of NaCl, KCl, MgCl 2 and CaCl 2 and a total salt content between 11.8 and 29.2 wt% were considered ( Table 11 ).As was the case with the MgCl 2 salt inhibitor, CSMGem could not provide predictions and was not considered in this test.MF RKSA and CSMHyd showed highly biased results with  of 41.5 and 25.8 MPa respectively.On the other hand, CPA based methods performed better with HF CPA keeping the lead with an || value of 10.9 MPa compared to the 18.4 MPa score of the MF CPA.

Detailed results
In this section a detailed analysis of the obtained results is attempted.For each system considered, the absolute pressure error (MPa) is correlated to the system temperature, the dissociation pressure itself, the inorganic component content and the number of data points.More specif-     ically, the temperature range is split in clusters each 2 K wide, and the average absolute pressure error is computed for each method tested only for the data points included in that cluster.The dissociation pressure and inorganic component concentration range, indicated by the green and the orange boxes respectively are also plotted against temperature.Each box extends from the lowest to the highest data value included in each corresponding cluster and the experimental data medians (50th percentile of pressure and component concentration) are indicated by markers plotted on each corresponding cluster box.The number of data points included in each cluster is also shown in the lower part of the Figure .This visualization, though full of information, gives the reader a clear idea of the number and variety of the tests performed in this work to evaluate the predictive accuracy of each package.Note that each pressure deviation data series has been slightly shifted horizontally within the temperature clusters to avoid overlapping and to enhance readability.

CO 2 --CH 4 mixtures
In the low temperature range of 274-282 K, all thermodynamic models predict hydrate dissociation pressure fairly accurately for temperatures below 282 K ( Fig. 1 ), exhibiting nearly identical behavior with an average absolute error of 0.35 MPa ( ∼50 psi).The average absolute error increases abruptly to 2 MPa for temperatures above 284 K (52°F) for all thermodynamic approaches tested.As both temperatures ranges share similar the CO 2 content range, the predictive accuracy of the models is mostly affected by the temperature increase, which also causes an increase in the dissociation pressure itself, rather than by the CO 2 concentration, verifying that system temperature is a dominant factor influencing the predictive performance.
HF CPA model shows slightly lower error increase with temperature, probably due to better tuned model parameters (i.e., Kihara parameters) for the CO 2 and/or CH 4 components for this type of fluids, or due to tuning using experimental hydrate dissociation data in a wider temperature range.On the other hand, CSMHyd exhibits poor performance compared to the other approaches tested, with the absolute error reaching values up to 2 MPa ( ∼300 psi).

CO 2 -N 2 -CH 4 mixtures
MF CPA, MF RKSA and HF CPA are shown to exhibit almost identical behavior over the entire temperature range with absolute errors less than 1.1 MPa, although pressure can be as high as 30 MPa ( Fig. 2 ).All three methods utilize the original VdW-P theory to model the fugacity of water in the hydrate phase but differ in the fugacity models (MF CPA/MF RKSA) used for the coexisting phases.
Although CSMGem exhibits fine behavior at temperatures below 285 K, larger deviations are observed for temperatures above that and pressures above 15 MPa.This can be attributed to the modified ver-sion of the VdW-P theory utilized by CSMGem and its limited ability to capture hydrate thermodynamics at such conditions.
Considering the effect of impurities concentration (N 2 -CO 2 ), the nearly same concentration range over all clusters (but the 285-287 K one) does not correlate to the observed deviations variance, thus demonstrating the dominant effect of dissociation temperature and the pressure values.

CO 2 -N 2 mixtures
All tested models except CSMHyd exhibit similar accuracy for the hydrate dissociation pressure of hydrate-forming CO 2 -N 2 mixtures in the 273-285 K range ( Fig. 3 ).CSMHyd, on the other hand, shows a linear increase in the deviation error with increasing temperature, regardless of the CO 2 content.This can be attributed to the fact that CSMHyd estimates dissociation pressure using experimentally calculated chemical potential difference rather than running flash calculations.Additionally, Langmuir constants are based on empirical relationships with fitted constant values which apply to specific temperature range and depend on the hydrate structure, the size of the cavity and the trapped gas molecule in hydrate cages (N 2 and CO 2 in this case).
All models except CSMHyd predicted hydrate dissociation pressure of CO 2 -N 2 mixtures in almost identical accuracy for systems with temperature below 285 K and high carbon dioxide concentrations (first four clusters with an average concentration of 80%).Pressure predictions are improved when the average CO 2 content drops down to 50%, indicating the importance of the CO 2 -N 2 content in the hydrate mixture.However, significant deviations are observed for the single data point available at high temperature which exhibits low CO 2 content and very high dissociation pressure value.The deviations should be attributed to the corresponding very high N 2 content and only MF RKSA, MF CPA and CSMGem exhibit errors less than 10% (5.8 MPa).

CH 4 -N 2 mixtures
All models share similar predictive performance over the entire temperature and pressure range (276-292 K and 5-25 MPa) with N 2 concentration and temperature effect being parameters that both affect predictions accuracy ( Fig. 4 ).Deviations exhibit significant correlation to the N 2 content due to their high values when the concentration is high, as well for the single data point in the 278-280 K cluster that exhibits very low N 2 content (10 wt%) and very accurate dissociation pressure prediction.Overall, temperature affects the predictive capability of the examined models as deviations exhibit an increasing trend with increasing temperature.As a result, the worst-case deviation of approximately 1.8 MPa is observed at the 286-288 K cluster when N 2 also exhibits high concentration.In total, all models still provide dissociation pressure predictions of high quality even when the pressure is as high as 20 MPa or more.

CH 4 -Saline aqueous phase mixtures (Pure salt: NaCl)
All tested models provide satisfactory results at low temperatures (272-282 K) and low NaCl concentration of < 8 wt% in the aqueous phase.However, when NaCl concentration increases (278-280 /16 wt%) all models show increased errors, with MF RKSA exhibiting highest deviation among all ( Fig. 5 ).The combination of moderate temperatures (282-290 K) to high NaCl content increases the dissociation pressure itself and introduces even more error to its predictions.At high temperatures, the performance is descent for all models with CSMGem and the CPA based models providing best results.
Interestingly the performance of the tested methods is very highly correlated over the full temperature and NaCl content with CSMGem providing remarkably accurate predictions for all temperatures and salt concentrations, even at dissociation pressures as high as 150 MPa, possibly due to the modification of VdW-P theory implemented.A more thorough comparison of the CPA-based models shows that HF CPA and MF CPA share similar deviation trend, while the latter provides more accurate predictions in the low to moderate temperature range (276-292 K) temperature range.The improved performance of MF CPA over HF CPA could be due to the addition of the virial coefficient and the born term to the Helmholtz residual energy of the system, which seems to be more representative of hydrate-forming systems with NaCl salt component.

CH 4 -Saline aqueous phase mixtures (Pure salt: KCl)
All tested models exhibit high accuracy in the low temperature range (275-283 K) with an average KCl concentration of 10-12 wt%, mainly due to the fact that the dissociation pressure to be predicted exhibits low values ( Fig. 6 ).When KCl concentration increases close to 30 wt%, CPA-based models suffer from an extraordinary error increase.For example, HF CPA and MF CPA show same significantly high deviations in As was the case with previous systems, MF RKSA and CSMHyd are not capable at all to catch the thermodynamics, while CSMGem provides best predictions for hydrate dissociation pressure among all approaches tested.

CH 4 -Saline aqueous phase mixtures (Pure salt: MgCl 2 )
The absolute pressure error (MPa) in the presence of an aqueous phase containing pure MgCl 2 depends clearly on the MgCl 2 salt concentration in the aqueous phase rather than on the temperature ( Fig. 7 ), as can be verified from the very strong correlation between the deviation error and the MgCl 2 content in the top two subplots.The deviations also correlate well to the dissociation pressure as an indirect result of the dissociation pressure correlation to the salt content.As mentioned in Section 4.1.7, CSMGem is not included as it cannot treat that specific salt.
HF CPA achieves the best by far predictive performance for both low ( < 10 wt%) and high MgCl 2 concentrations ( > 10 wt%) data points in all temperature ranges, as far as both the absolute and the relative pressure error are concerned.Although MF CPA exhibited very similar accuracy for uninhibited systems, it now exhibits significantly larger deviations in the present system (3-6 times higher deviation than HF CPA).This needs to be attributed to the poor performance of the virial corrections and Born term to the residual Helmholtz energy that describe this specific salt effect.On the other hand, HF CPA performs better thanks to the properly tuned Debye-Huckel term, which accounts for the effect of salt inhibition on hydrate dissociation conditions.

CH 4 -Saline aqueous phase mixtures (Pure salt: CaCl 2 )
When it comes CaCl 2 inhibited systems, the need to interpret results according to their temperature, discussed in Section 4.1.8, becomes now clear ( Fig. 8 ).Clearly, all models exhibit very poor predictive performance when extrapolating at such low temperatures (257-265 K) and CaCl 2 salt concentrations averaging to 30 wt%.On the other hand, all models perform significantly better at moderate and high temperatures although those clusters share similarly high salt content and dissociation pressure values.
All tested models show very good results at moderate temperatures, regardless of the salt concentration.MF RKSA exhibits slightly higher, still affordable, deviations and that is attributed to the NaCl pseudo component method used to describe the salt content in the aqueous phase.At high temperatures, above 283 K, the CPA based models and CSMGem provide very accurate predictions even at higher CaCl 2 concentration (23 wt%).MF RKSA still suffers from the same drawbacks which now lead to extremely high deviations of 50 + MPa while CSMHyd provides large but reasonable errors.

CH 4 -Saline aqueous phase containing salt mixtures
The absolute pressure error for systems inhibited by salt mixtures (e.g., MgCl 2 + NaCl) is shown in Fig. 9 .Results indicate that for low temperatures (below 274 K) and salt concentration of 29 wt%, all models except HF CPA deviate severely by producing errors of more than 50 MPa, whereas the latter shows very good performance with errors of less than 10 MPa against dissociation pressures of more than 100 MPa.As expected, MF RKSA keeps failing to accurately predict the dissociation pressure even at moderate to high temperatures while CSMHyd exhibits better performance which, however, leads to very large deviations at high temperatures.This could be attributed to the switch of the Langmuir constant evaluation method utilized from the two parameter expression for temperatures below 300 K to the Kihara potential based one with innacurate model parameters for temperatures above that.On the other hand, the CPA based models perform much more accurately in the moderate and high temperature range with HF CPA ranking first, despite the high concentration of the salt mixture in the aqueous phase (up to 25%) and the high dissociation pressure values.

Comparison
Arriving to concrete conclusions on the software packages prediction capability is not an easy task as distinct strong points and weaknesses have been identified for each individual system in the previous Section.To attempt that, two plots were generated.Firstly, the average behavior of each software package was further evaluated by comparative plots showing the  ±  values.Fig. 10 demonstrates the errors obtained for all uninhibited systems and it is clear that the software packages tested exhibit significantly variant performance.CH 4 -CO 2 mixtures are mod-eled accurately with small only deviations which can be rendered as affordable for the handling of flow assurance issues in pipeline networks.Slightly enhanced, still affordable errors are obtained for the CH 4 -CO 2 -N 2 and the CH 4 -N 2 systems as the maximum average absolute error does not exceed 1 MPa (145 psi).However, this is not the case with CO 2 -N 2 mixtures where all approaches exhibit weak performance with errors equal to or above 1.7 MPa.Interestingly, all software packages examined exhibit similar performance for each single system considered (similar colors in each row), thus rendering the prediction capability of the available modeling techniques as rather uniform.
A different picture is obtained when considering inhibited systems where results similarity is observed over packages (similar colors in each column) rather than over data systems ( Fig. 11 ).Errors are not only significantly enhanced but they also exhibit great variance among the systems examined.CSMGem provides predictions of superior accuracy compared to the rest of the packages for systems it is capable handling   of (it cannot handle MgCl 2 ).On the other hand, MF RKSA provides notoriously large deviations, especially for systems containing MgCl 2 and CaCl 2 thus demonstrating the weakness of the pseudo-salt NaClequivalent component approach.Although CSMHyd exhibits better performance than MF RKSA it still provides highly deviating predictions for most of the inhibited systems examined.Both CPA based models exhibit descent performance over all systems thus providing a safe dissociation pressure prediction solution for inhibited systems.
To further facilitate the performance comparison, a set of parity plots has been added to the supplementary material demonstrating the exact performance of each system.Interestingly, specific data points exhibit same behavior when tested on all five packages.As an example, consider the bunch of five measurements on uninhibited systems containing N 2 and CO 2 at pressures between 2.4 and 3.2 MPa (Fig. S3) performing so similarly between all five packages, that it might need to attribute them to experimental rather than modeling errors.The weakness of MF RKSA when it comes to strongly inhibited systems is apparent in Figs.S8 and  S9.
An overall evaluation can be attempted by means of the spider plots in Fig. 12 .For classic, uninhibited systems, all packages provide very similar, quite unbiased results with average absolute errors barely exceeding 2 MPa, thus leading to no winners.When it comes to inhibited systems, MF RKSA and CSM Hyd exhibit severe errors whereas HF CPA and MF CPA perform much better with HF CPA taking the lead due to its inferior and better balanced average error.CSM Gem exhibits optimal, best balanced performance but it suffers from its inability to handle systems such as those inhibited with MgCl 2 .
Although pressure errors for fixed temperature have been considered in this work thanks to their direct relation to design decisions, when huge deviations are observed they may correspond to minor only errors in the estimation of the related dissociation temperature for fixed pressure.This is the case when the dissociation conditions considered lie very close to the maximum hydrate forming temperature where the hydrates curve shape usually exhibits a very steep shape, nearly vertical shape in the p-T plot.As a result, software users are encouraged to run a very simple sensitivity analysis when requesting dissociation pressure estimates, by varying the temperature and monitoring the change in the predicted pressure.A huge sensitivity (i.e., derivative) value would indicate that the examined conditions are very close to the steep part of the dissociation curve where predictions are unreliable, hence they should be avoided in design tasks.

Conclusions
In this work, five commercial software packages (HF CPA, MF CPA, MF RKSA, CSMGem and CSMHyd) for the prediction of hydrate dissociation conditions were presented, investigated, and compared against a large experimental database of inhibited and uninhibited 400 data points.The results of this study show that all five models provide reliable, balanced predictions of dissociation pressure for uninhibited flow assurance systems within an error range which is considered as acceptable to ensure successful and economical hydrate-free gas flow.HF CPA provided slightly more accurate predictions, while all models accuracy was shown to be sensitive to temperature rather than to the concentra-tion of impurities (CH 4 -CO 2 and CH 4 -CO 2 -N 2 systems).For the technological applications collection (CO 2 -N 2 and CH 4 -N 2 systems), all models exhibited more severe sensitivity of the prediction accuracy to N 2 concentration rather than to the system temperature.It was noticeable that all models, except CSMHyd, showed almost identical performance and the predictions were considered sufficiently reliable for most technological applications.
For inhibited systems with monovalent cations (NaCl, KCl), CSMGem was found to exhibit remarkably balanced accuracy for all temperatures and salt concentrations compared to other models tested.CPA-based models (MF CPA, HF CPA) are ranked as the second-best predictors, although they cannot cope with in high salt concentrations ( > 25 wt%).MF RKSA and CSMHyd exhibited poor predictive performance and were unable to determine dissociation pressure in descent accuracy.For systems with divalent cations (MgCl 2 , CaCl 2 ), HF CPA was found to provide remarkably accurate predictions for both salts while MF CPA performed accurately only for CaCl 2 .CSMGem produced dissociation pressures for CaCl 2 with moderate accuracy and showed no capability of handling MgCl 2 , while MF RKSA and CSMHyd performed poorly, with MF RKSA achieving poorest performance.
For the gallery of methane hydrate-forming systems with salt mixtures, which is a typical case for dry gas production streams, CPA-based methods demonstrated best performance, with HF CPA providing the most accurate prediction.In contrast, MF RKSA and CSMHyd provided very poor predictive performance and were unable to predict accurately dissociation pressure values for flow assurance or other potential hydrate applications (e.g., desalination).
Overall, all software packages provided similar performance, sufficiently reliable to obtain dissociation pressures for the uninhibited systems.However, when considering results on inhibited systems, each package exhibits its own performance profile with individual strengths and weaknesses.It could be said that HF CPA ranks as the most accurate tool in predicting dissociation pressure compared to other models tested due to its balanced behavior over all types of systems followed closely by MF CPA and by other models in the order CSMGem, CSMHyd and MF RKSA.
From a thermodynamic point of view, our results indicate that the safer solution currently available to go with hydrate dissociation conditions prediction is the combination of the Van der Waals-Platteeuw theory with a regular fluid phase model, such as the SRK, and an advanced fluid model which also considers cross-association hydrogen bonding and electrolytes.Nevertheless, this approach does not constitute a predictive model implying that parameters tuning, such as the Kihara potential terms, needs to be updated in the arrival of new experimental data.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 9 .
Fig. 9. Detailed error distribution of the CH 4salt mixture system.

Fig. 12 .
Fig. 12. Overall MAE and ME per system and per software package. .100103

Table 3
Performance of the CO 2 -CH 4 system.

Table 4
Performance of the CO 2 -N 2 -CH 4 system.

Table 5
Performance of the CO 2 -N 2 system.

Table 6
Performance of the CH 4 -N 2 system.

Table 7
Performance of the NaCl-salt system.

Table 8
Performance of the KCl-salt system.

Table 9
Performance of the MgCl 2 -salt system.
* CSMGem does not support systems containing MgCl 2 .

Table 10
Performance of the CaCl 2 -salt system.

Table 11
Performance of salt mixture systems.