Applying SAFT-type models for the anomalous properties of water: Successes and challenges

Water is the most important, but at the same time the most anomalous substance on earth. Due to its anomalies, it is a challenge even for advanced equations of state like PC-SAFT, CPA and SAFT-VR Mie to predict water’s anomalous properties. It is believed by some researchers that water’s anomalous properties are the result of structural fluctuations. The recently developed DAPT model is consistent with this idea and it is actually able to predict a considerable amount of water’s anomalies, like the density maximum and speed of sound maximum, which is promising, but we show in this work that it produces also some unrealistic results like minima and maxima for isochoric heat capacity. Also DAPT’s accuracy falls at high temperatures/high pressure regions where water behaves like a simple fluid. This reveals further the complexity of the problem since a successful model should treat water as an anomalous and as a simple fluid when necessary. Moreover, all of the association models fail to reproduce some of water’s structural properties and especially the tetrahedral fraction. It may be that a more accurate structural representation of water and the inclusion of steric hindrance and cooperativity effects could provide a more realistic model for water.


Introduction
Water is the most important liquid and at the same time the most anomalous one [1].Most of water's thermodynamic properties show extrema as temperature and/or pressure change.A well-known example is the density maximum at 4 • C at atmospheric pressure.Due to its anomalous thermodynamic behavior, it is a challenge for any theoretical thermodynamic model to accurately describe its properties [2].Some extremely noteworthy thermodynamic models are those that have been derived from the statistical associating fluid theory (SAFT) [3,4], which has a more realistic approach on modeling intermolecular interactions.As a result, it is more advanced compared to classical equations of state (EOS), but it is not the only advanced one.Another great advantage is that it considers hydrogen bonding (association) between molecules explicitly.The contribution of association forces is based on the Wertheim's perturbation theory [5][6][7][8].
Throughout the years many different EOS have been developed that are based on SAFT.SAFT and Wertheim's formalism rely on the structure and thermodynamic properties of a reference system.The original SAFT considers a hard-sphere reference system [3,4].Some SAFT variants are based on a different reference system.For instance, the Perturbed-Chain SAFT (PC-SAFT) [9] is based on a hard-chain system and soft-SAFT [10] is based on a system of Lennard-Jones chains.
Other EOS have completely different expressions for the physical forces and they have adopted an association term extremely similar to that of SAFT, like the Cubic-Plus-Association (CPA) EOS [11].There are also some variants to SAFT and Wertheim's perturbation theory that incorporate more phenomena and interactions like polar forces [12,13] and cooperativity [14,15].
These were a few examples of the various SAFT-type models.The list of all the different SAFT-type models is quite long with many reported successes of the models.For example, PC-SAFT and CPA can be very accurate at describing the solubility of water in hydrocarbons [16][17][18][19] and many SAFT-type models are accurate at describing vapor-liquid equilibria of alcohols with hydrocarbons [17,20].However, in the case of pure liquid water none of these aforementioned models is able to capture most of its anomalous thermodynamic properties [2].It seems that all of these models do not account for the origin of water's anomalous properties (whatever this may be).However, the true origin of water's anomalies is still a matter of debate in the scientific community [1,21].
Some researchers argue that the origin of water's anomalies is structural fluctuations.In specific, they argue that there are fluctuations between two different local structural environments in water [1,22,23].Marshall has also stated that the water's anomalous behavior is because https://doi.org/10.1016/j.fluid.2022.113617Received 8 July 2022; Received in revised form 19 September 2022; Accepted 21 September 2022 of its transition into tetrahedral symmetry [24].In a series of publications [25,26], Marshall has introduced changes to the structure of the reference system so that it would resemble more accurately water's structure.More specifically, tetrahedral coordination is forced to fully hydrogen bonded molecules.Interestingly enough, with these changes the doubly associated reference perturbation theory (DAPT) was able to capture some of water's anomalies, like the density maximum and the heat capacity minimum with respect to temperature.
This two-state hypothesis has been supported by molecular simulation models [22,27] and has been implemented in other thermodynamic models (mostly by the group of Anisimov) [28,29].These thermodynamic models are very successful at capturing water's thermodynamic anomalies.These models however are semi-empirical with many adjustable parameters and they have been fitted to many of water's thermodynamic properties.When they are applied at high temperatures the models produce very unrealistic results [30].It has been discussed that all these thermodynamic models that are consistent with this hypothesis are quite successful at describing water's anomalies in contrast to models that do not follow it [27].
As it has been discussed by Tsochantaris et al. [2], water's thermodynamic behavior is largely dependent on pressure.For example, density no longer shows maxima at high pressures (>200 MPa).Thus, it is important to evaluate a model's performance over a wide range of conditions.As a continuation of a previous study [2], we are going to apply DAPT for water and compare it with other advanced equations of state for a wide range of conditions and for many thermodynamic properties.We are also considering some important structural properties of water, like the fraction of molecules that are tetrahedrally coordinated.The other models that we are going to consider are CPA, PC-SAFT and SAFT-VR Mie (variable range SAFT for Mie fluids).
Nezbeda et al. [31] have very recently published a comprehensive study where they have applied multiple SAFT-type models for water.They have considered PC-SAFT, CPA and SAFT-VR Mie as in this work, but they also considered more models like polar PC-SAFT with different polar terms [32,33], critical point PC-SAFT [34] and others.In addition, they considered the association dependent PC-SAFT (ADPC-SAFT) [24], which proved to be the most accurate model according to their evaluation.The key change in ADPC-SAFT is that the segment diameter is dependent on the amount of fully hydrogen bonded molecules.While ADPC-SAFT shows great accuracy for many properties, it is not successful at capturing some of water's anomalies like the density maximum.It is, however, able to produce heat capacity minima.Our study focuses especially on water's anomalies along with some structural properties that are considered to be connected with water's anomalies [1,27].For this reason, we have considered DAPT which is able to produce various anomalies (more than any of the other aforementioned models).
In the next section of the manuscript (Models and Database), the thermodynamic models used in this work are presented along with expressions for thermodynamic and structural properties relevant to this work.In the same section, an overview of the experimental data for the thermodynamic properties is shown.In the section ''Results and Discussion'', results from the models are shown and compared with data from the literature.Finally, the conclusions of this work are presented.

PC-SAFT
The PC-SAFT equation of state was developed by Gross and Sadowski [9] and is a modified version of the original SAFT model where the reference system is a hard-chain system.In this study, the simplified PC-SAFT proposed by von Solms et al. [35] is used.For pure components the simplified version has been reported to produce the same results as the original [9].For the application of PC-SAFT, five adjustable parameters need to be specified.These are usually fitted to vapor pressure and saturated liquid density data.The reduced residual Helmholtz free energy in PC-SAFT can be expressed from the following equation [16]: where  ℎ is the Helmholtz free energy contribution of the hard sphere interactions,  ℎ is the contribution from chain formation,   is the contribution from the dispersive forces between whole chains of hard sphere segments,   is the contribution from association forces (hydrogen bonds) between segments and  is the chain length (segment number) of the molecule.All the energy contributions are dimensionless.
The hard sphere term in the simplified version is expressed as: where  is the packing fraction and it can be expressed as: where  is the number density,  is the chain length (segment number) of the molecule and  is the temperature dependent diameter.The chain term  ℎ for a pure fluid can be expressed as: where  ℎ is the radial pair distribution function for hard sphere segments of the component.The dispersion forces are evaluated using second order perturbation theory [36,37]: where  1 is the first order term and  2 is the second order term.These terms for a single component can be expressed as [38]: where  is the segment energy parameter,   is the Boltzmann constant,  the temperature of the system and  ℎ is the hard chain contribution to the pressure  .The terms  1 and  2 are the reference systems integrals with: where  = ∕ is the reduced distance and  is the range of the squarewell interactions.In order to estimate the integrals  1 and  2 , Gross and Sadowski [9] used the polynomial terms: where   () and   () are empirical terms that are dependent on the chain length .The universal coefficients of these terms can be found in the original publication [9].The association term   for a single component can be expressed as: where   is the fraction of the molecules not bonded at site  and  is the number of association sites on the molecule.It is common to assume that all hydrogen bonds have the same association strength and as a result in some cases, like the 4C scheme, every site has the same free site fraction   .
The site fraction   can be expressed by: (10) where   is the association strength between two sites.More extensive information on the PC-SAFT equation of state and the simplified PC-SAFT is available in the original literature [9,35].

CPA
Kontogeorgis et al. [11] developed the Cubic-Plus-Association (CPA) EOS by combining the classical Soave-Redlich-Kwong (SRK) EOS with an association term very similar to the association term of the SAFT EOS.As a result the expression for the physical forces are much simpler compared to those of the SAFT EOS.For the application of CPA, five adjustable parameters need to be specified.These parameters are typically fitted to vapor pressure and saturated liquid density data.The reduced residual Helmholtz free energy of the system can be expressed as [39]: where   is contribution of the SRK EOS to the reduced residual Helmholtz free energy and   is the contribution of association interactions (hydrogen bonds).
The  SRK can be expressed as [39]: where  is the co-volume parameter,  is the molar volume, ( ) is the temperature dependent energy parameter.The contribution of association forces   is very similar to the association term of PC-SAFT (Eq.( 9)).The main difference resides in the calculation for the association strength   since there is a different expression for the radial distribution function g.More information on the CPA is available in the original literature [11].

SAFT-VR Mie
Lafitte et al. [40] developed SAFT-VR Mie where the reference (monomer) system consists of sphere segments that interact with a Mie potential (generalized Lennard-Jones potential) [41].In SAFT and PC-SAFT, the dispersion forces are expressed via a simple square-well potential.As a result, SAFT-VR Mie has much more complex terms for the dispersion forces, which will not be reviewed in detail here.For the application of SAFT-VR Mie, seven adjustable parameters are required.The reduced residual Helmholtz free energy can be expressed as: where   is the Helmholtz free energy contribution of the monomer (reference) fluid.The contributions  ℎ and   have the same meaning as the ones in PC-SAFT.The monomer contribution   is a sum of the hard sphere contribution and dispersion forces contribution (  =  ℎ +   ).The monomer contribution for a single component can be expressed as: where  is the chain length (segment number) of the component,  ℎ is the hard sphere contribution (same as the original SAFT),  1 ,  2 ,  3 are the first, second and third order perturbation terms respectively and  = 1∕(   ).
For the association term, we have used the simplified LJ-based kernel which was introduced by Dufal et al. [42].More details about the expressions and the model in general can be found in the original publications [40,42].

DAPT
In most SAFT variants, the structure of the reference fluid is considered to be unaffected by hydrogen bonds.While this might be accurate for many substances, it is certainly not accurate for the case of water.Because of hydrogen bond formation, some water molecules are tetrahedrally coordinated.In DAPT [26], the structure of the reference fluid is affected by hydrogen bonds and tetrahedral symmetry is enforced at molecules that are fully hydrogen bonded.It should be noted though that DAPT was developed only for water and was not extended to mixtures.The reduced residual Helmholtz free energy of the s ystem is: For the  ℎ , the Carnahan-Sterling equation was used [43].The dispersion forces are expressed via a square well potential of range .The   is expressed with a first-order perturbation theory [36]: where  is the number density of molecules ( = ∕ ) and   () is the integral of the reference system correlation function over the range  of the square well potential.This dispersion term is essentially the same as the first-order perturbation term of the dispersion contribution in the Associated Perturbation Theory (APT) [25].As commented by Marshall [26], a second-order perturbation term does not improve the methodology and therefore it has been omitted.The   () is expressed as : where   () is the correlation function of the reference fluid.In the original SAFT [4], the reference system fluid is a hard-sphere fluid with   =  ℎ , while in PC-SAFT [9], the reference system is a hard-chain fluid with   =  ℎ (Eq.( 7)).In DAPT, the structure of the reference system is different and it is considered to be affected by hydrogen-bonds (association), as previously mentioned.The   in DAPT is expressed as: where  4 is the fraction of water molecules bonded 4 times and  ℎ (  ) is the integral of the correlation function of the hard-sphere system over the critical radius   .The radius   is the range of the association forces (i.e.H-bonds cannot be formed at distances larger than   ).Eq. ( 19) is derived by enforcing tetrahedral coordination at molecules that are hydrogen bonded 4 times (i.e. by forcing the coordination number to be equal to 4).Through this change, there are essentially two structures in the reference system.The molecules that are bonded 4 times follow the tetrahedral coordination, while the rest have the structure of a hardsphere fluid (like in the original SAFT EOS).If there are no molecules that are fully hydrogen bonded ( 4 = 0), then the integral   is the same as the hard-sphere integral  ℎ .It should be noted that this change to the integral affects the dispersion term as well as the association term.More information about the derivation of Eq. ( 19) can be found in the original publications [25,26].
The association term of DAPT by assuming the equivalence of sites is expressed as: where   is fraction of the unbonded molecules and  (0) is the graph sum that contains all the information of the hydrogen bonding interactions.Note that this term is significantly different from the one in Eq. ( 9).Because the integral of the correlation function is dependent on  4 , the association term needs to be derived from the more general Wertheim term, as demonstrated by Marshall [26].
The graph sum  (0) is expressed as : where   is the probability of two association sites are oriented correctly to connect (also referred as association volume) and   is the magnitude of the Mayer function of the site-site association interactions.Both of these properties are connected to the adjustable parameters.
It should also be noted that to estimate the unbonded fraction   one will have to solve the following equation: where   is: where   − =   is the density of all molecules that are not bonded at site A.

Thermodynamic and structural properties
In this study, many thermodynamic properties will be considered that are connected to water's anomalous properties.The expressions for the isothermal compressibility   , isobaric thermal expansivity   , residual isobaric heat capacity    , residual isochoric heat capacity    and speed of sound  are as follows: It should be noted that for some calculations the ideal gas contributions of the heat capacities is needed.Values for these ideal gas contributions are obtained from the DIPPR database [44].
Besides thermodynamic properties, the SAFT-type models can also evaluate some structural properties as well.An important structural property for associating compounds like water is the k-times bonded fractions.In the SAFT framework, it is common to assume that each site in a four sited molecule has the same probability to form a hydrogen bond (H-bond).At the same time, each hydrogen bond formation is independent from the rest (association probability/volume and energy are the same).Due to these characteristics, the probability of forming  hydrogen bonds can be estimated by using a binomial distribution [45].The binomial experiment has 3 key characteristics: (a) A fixed number of ''attempts'', which in this case is equivalent to the number of sites that can form a bond, (b) there are only two possible outcomes with ''success'' and ''failure'', which in this case are equivalent to ''form a Hbond'' and ''not form a H-bond'' and (c) the ''attempts'' are independent and have equal probability at occurring and this is obviously true with the equivalence of sites in the 4C scheme where each site has the same probability to form a H-bond.Thus, the fraction of k-times bonded molecules can be expressed as [45]: where   is the fraction of molecules not bonded to site A and   is the number of sites.

Model implementation
The implementation of DAPT was accomplished with the use of Clapeyron.jl[47].Clapeyron.jl is an open-source software that contains a library of many published thermodynamic models and especially SAFT-type models and at the same time a useful framework that for the implementation and development of many thermodynamic models.The library includes already PC-SAFT, CPA and SAFT-VR Mie, which have also been used in this study.One of the key advantages of Clapeyron.jl is that it relies on automatic differentiation for the calculation of derived properties.As a result the user has to implement only the terms of the Helmholtz free energy for the calculation of all the thermodynamic properties.For the automatic differentiation, Clapeyron.jlrelies on the package ForwardDiff.jl[48].

Database
To evaluate the performance of the various thermodynamic models, experimental data have been collected for various thermodynamic properties.For some properties and for some conditions, experimental data are scarce.For this reason, values from the NIST Chemistry Webbook [49] were used as well.These values are generated from the IAPWS-95 formulation for ordinary water [50] which is a highly accurate multi-parameter model for water.An overview of the data and the values from NIST is shown in Table 1.Of all the different data, it is important to comment on the values from Troncoso [51].The heat capacity values of Troncoso are significantly different from the ones that are estimated from NIST as shown in Fig. 1.At high pressures, the reported qualitative behavior is completely different.The IAPWS-95 formulation [50] suggest that the heat capacity shows an increase even leading to maxima.A similar behavior has been reported by Lin and Trusler [52].However, the data by Troncoso [51] suggest a different trend with no maxima at these conditions.It is worth noting that the IAPWS-95 formulation was not fitted to any   data at low temperatures (275-300 K) and high pressures (200-400 MPa).Troncoso estimated the   by using a calorimeter.For this study, we are going to consider both datasets for the evaluation.

Deviations
To estimate the quantitative accuracy of the various models, we shall use the absolute average deviation (%AAD).The %AAD of an estimated property  is expressed as: where  , and  , are any thermodynamic or structural property at the experimental point  calculated from the models and estimated from experiments respectively and   is the total number of experimental points.Note: To calculate the residual isobaric and isochoric heat capacity from the data [51,52] and the values from NIST [49], the ideal gas contribution from the DIPPR database[44] was used.Note:   and   are saturated pressure and saturated liquid density respectively. is the co-volume parameter,   is association energy,   is the association volume parameter and  is the ideal gas constant.Finally,  =  0 ∕() and  1 and  0 are constant parameters related to the temperature dependent ( ).Note:   and   are saturated pressure and saturated liquid density respectively. is the chain length,  is the temperature independent segment diameter,   ,   are the exponents characterizing the repulsive and attractive dispersion interactions respectively,  is the segment energy parameter,   is association energy and   is the association volume.Note:     are saturated pressure and saturated liquid density respectively. is the temperature independent segment diameter,  is the range of the square-well interactions,  is the segment energy parameter,   is the association energy,   is the critical angle and   is the critical radius of association interactions.

Thermodynamic properties
The application of the thermodynamic models requires the specification of a few adjustable parameters.The parameters are typically adjusted to vapor pressure and saturated liquid density data.For PC-SAFT and CPA, there are plenty of parameter sets available in literature for water and there are some that consider water molecules to form up to 2 bonds (2B scheme), up to 3 bonds (3B scheme) and up to 4 bonds (4C scheme).In this work, we are going to consider parameter sets that only follow the 4C scheme, which is more realistic since it is generally known that water can form up to 4 hydrogen bonds.Also, this selection will be important to estimate the fraction of molecules that are 4-times hydrogen bonded.The parameter sets of each model are shown in Tables 2-5.All of them have been fitted to saturation pressure and saturated liquid density data though at different temperature regions.Few of the investigated parameter sets considered phase equilibria of water containing mixtures [16,61].All of the parameter sets were applied for the thermodynamic properties that are shown in Table 1.The comparison results for the various thermodynamic properties are presented in Figs.2-6.
From Fig. 2, it is clear that of all the models investigated, only DAPT is able to predict the density maximum and in fact it shows exceptional accuracy for water at 0.1 MPa.At higher pressures, the accuracy changes however.Water's thermodynamic behavior is heavily affected by pressure to the point where at high pressures (≥200 MPa) there are no density maxima.At these conditions, other SAFT-type models actually perform better.The comparison in Fig. 2 shows the great complexity of water's thermodynamic behavior since a model should not only be able to capture the density maximum, but it should also revert back to normal fluid behavior at high pressures.Due to this shift in the thermodynamic behavior, the set P1 (PC-SAFT) actually has the highest accuracy at high pressures.
In Fig. 3, the isothermal compressibility and the isobaric expansivity of liquid water is displayed.Another anomaly of water is the isothermal compressibility minimum at approximately 315 K.It has been noted that this minimum does not shift with pressure [62].With the exception  [49,50] and symbols are experimental data by Troncoso [51].Fig. 2. Density isobars of liquid water at different conditions.Lines are predicted results from the models and symbols are data and results from NIST [49,50] (Table 1).
of DAPT, all other SAFT-type models do not capture any minima.DAPT is able to capture minima but at higher temperatures (approximately 370 K).The isobaric expansivity of water does not show any extrema, but it is directly connected to the density maxima where   = 0.Only DAPT is able to predict values with   ≤ 0 and it even shows very great accuracy for 0.1 MPa (Fig. 3).However, the accuracy of DAPT significantly decreases at high temperatures (especially at 400 MPa).This is obviously connected to the lower accuracy of density estimation at 400 MPa (Fig. 2) where the expansion rate of water is overestimated.On the other hand, the other models perform much better at high pressures and at all pressure points the isobaric expansivity isobars of these models (CPA, PC-SAFT, SAFT-VR Mie) have similar curvature.
In Fig. 4 estimations for the residual isobaric heat capacity are shown.As it was discussed earlier, in this work two contradicting datasets have been used and the model predictions are compared against both sets.The key difference in both datasets is that the heat capacity values are significantly different at high pressures as it was discussed before.For the estimation of the experimental residual heat capacity, we used total isobaric heat capacity values [49,51,52] and ideal gas contribution estimates from the DIPPR database [44].From Fig. 4, it is clear that only DAPT is able to predict some extrema.At 0.1 MPa, DAPT predicts a minimum at around 325 K, while the true one is at around 313 K. Also, the predicted minimum is more pronounced and the heat capacity values are underestimated.The predicted values of the rest of the models are also underestimated at 0.1 MPa.At higher pressures the accuracy of CPA, PC-SAFT and SAFT-VR Mie improves significantly (for both datasets).In fact, CPA and SAFT-VR Mie show excellent accuracy at high temperatures and high pressures (Fig. 4a).However, at these conditions DAPT predictions deviate more from the data and they even exhibit maxima.This unusual behavior was also noted by Marshall [24].None of the models show maxima at low temperatures and high pressures similar to the values by NIST.Similar conclusions are drawn, when the model predictions are compared against the data of Troncoso (Fig. 4b).All the model predictions underestimate the residual isobaric heat capacity values, but qualitatively they are more similar to the data [51] since these data suggest that there are no maxima.
Fig. 5 shows results for residual isochoric heat capacity and speed of sound.For the estimation of the experimental residual heat capacity, we used total isobaric heat capacity values [49,51,52] and ideal gas contribution estimates from the DIPPR database [44].Interestingly enough, the residual isochoric heat capacity does not show any extrema when plotted against temperature at 0.1 MPa.There are maxima though at 200 MPa and 400 MPa.In general, the model predictions deviate significantly from the actual values.The model predictions that show maxima are those of CPA, SAFT-VR Mie and DAPT.The predicted maxima though appear at much higher temperatures.Additionally, DAPT predicts minima as well, which is unrealistic for the    of water.The speed of sound of water shows maxima at all pressure points.Again only DAPT is able to predict anomalous behavior, where calculated speed of sound increases with an increased temperature and at high pressures (≥100 MPa) it predicts maxima.These maxima, however, are at higher temperatures compared to the actual ones.The predicted results of PC-SAFT, CPA and SAFT-VR Mie show similar trends where the values continuously decrease when the temperature increases.As a result, their calculated values deviate significantly at 0.1 MPa.At high pressures and temperatures, the accuracy of the models is improved.
In Fig. 6, the enthalpy and internal energy of vaporization is displayed.In general, the model predictions do not deviate significantly from the actual values.The most accurate models for these properties seem to be CPA and SAFT-VR Mie, while DAPT shows the largest deviations especially at high temperatures.
In order to quantify the accuracy of the various models, we use the %AAD.The %AAD for each of these properties was estimated using Eq. ( 31) and the results are shown in Table 6.As discussed earlier, two different data sets for the isobaric heat capacity have been considered and comparisons against both are presented.From Table 6, it is clear that no model (or parameter set) shows excellent accuracy for all properties.The parameter set P1 (PC-SAFT) performs best for density (quantitatively) at both unsaturated and saturated conditions.This is quite interesting considering the fact that PC-SAFT is unable to predict density anomalies.Meanwhile DAPT, which is known to be able to predict density anomalies, performs quite well quantitatively for density where the %AAD values are very close to those of set P1 (PC-SAFT), but it is not the best model.The set P1 (PC-SAFT) is also quantitatively the best for the isothermal compressibility   .The isothermal compressibility is directly related to the pressure derivative of density.It is likely that the set P1 (PC-SAFT) performs the best for these properties because it was fitted to a wider temperature range (275-640 K) compared to the rest of the parameter sets.The models show the greatest deviations for the isobaric expansivity.This is to be expected since it is a property directly related to the density anomalies.[49,50] (Table 1).

Fig. 4.
Residual isobaric heat capacity isobars of liquid water at different conditions.(a) Symbols are data [52] and results from NIST [49,50] (b) Symbols are data from Troncoso [51].Lines in all subfigures are predicted results from the models.For the ideal gas contribution, the DIPPR database [44] was used.
Only two models (DAPT and SAFT-VR Mie) show %AAD values below 100%.Even though DAPT is able to capture the anomalous expansivity, SAFT-VR Mie performs much better at high pressures where water behaves like a simple fluid for density.For the isochoric heat capacity, most model results deviate more than 10% with the exception of SAFT-VR Mie (9%).For the isobaric heat capacity (NIST values), most models also show deviations larger than 10%, where only CPA and SAFT-VR Mie show %AAD values smaller than that with 7% and 6% respectively.Deviations for all models are larger when model results are compared to data from Troncoso [51] because these data cover a more narrow range (275-310 K) and in that range of low temperatures the heat capacity tends to be anomalous.For   ,   and   , all models show satisfactory results with CPA showing the lowest deviations for these properties.Finally, the average %AAD for all the properties was estimated by averaging over the %AAD values for each property.For     there are two conflicting datasets and for this reason the calculation is performed twice where each time a different dataset for    is used.It turns out that SAFT-VR Mie provides the lowest average of deviations and the second best is DAPT.The average deviations of CPA (22%) are very close to those of DAPT (21%).We should state that there are other  [49,50] (Table 1).Fig. 6.Heat of vaporization (left) and Internal energy of vaporization (right) of liquid water.Lines are predicted results from the models and symbols results from NIST [49,50] (Table 1).models that could produce more accurate results like ADPC-SAFT [24] as discussed by Nezbeda et al. [31].One should consider, however, that during the parameterization of ADPC-SAFT more properties were included in parameter estimation procedure, e.g. the isobaric heat capacity.So the higher accuracy may also be, at least partially, due to the fitting procedure.

Structural properties
We are considering now some structural properties.It is a common practice for authors to evaluate association models for the fraction of free site   .For the case of water, one should also consider the fraction of molecules that are tetrahedrally coordinated   .For all of the association models, we are going to consider that the amount of tetrahedrally coordinated molecules are the same as the amount of molecules that are fully bonded  4 =   .In reality this approximation might not be entirely true since forming four hydrogen bonds is not on its own enough for tetrahedral coordination.However, the amount of molecules that have the correct arrangements cannot be estimated by our models since they do not include such high structural detail.For this property we have considered estimates from computational studies as well as experimental ones.Pathak et al. [67], and Russo and Tanaka [66] based their estimations on molecular simulations and they used the same method for identifying the tetrahedrally coordinated local structures.The only difference was in the model where Pathak et al. [67] used ST2 model and Russo and Tanaka used TIP5P [66].This suggests that the two models could have some differences at describing water's structure (i.e.different radial distribution function).Muthachikavil et al. [68] used molecular simulations of the iAMOEBA water model to identify local structures.They used a method similar Notes: The lowest %AAD values at each row are in bold.In the last two rows the average %AAD of all the properties is estimated.Since there are two conflicting datasets for the    , the average %AAD for all the properties is estimated two times by including a different    dataset.to the one described by Russo and Tanaka [66] but the criterion for the two-states is different since it is defined by the average angles of oxygen atoms with all its neighbors.Mallamace et al.Our results are shown in Fig. 7.The estimates of experimental studies and simulations show some differences quantitatively.Qualitatively though they tend to have the same sigmoid trend.Results from the thermodynamic models on the other hand deviate significantly from the other estimates quantitatively as well as qualitatively.Quantitatively most of the models largely overestimate the amount of tetrahedrally coordinated molecules.DAPT obviously shows the lowest deviations for all the values that were considered.Qualitatively, all the thermodynamic models show linear trends for this property.Other properties that were considered at 1 bar are the monomer fraction  0 and the fraction of partially hydrogen bonded molecules    =  1 +  2 +  3 .Results for both of these properties are shown in Figs. 8 and 9. Again, DAPT estimations are the ones closer to the experimental data.They still deviate significantly from the reported Fig. 10.Fraction of k-times bonded molecules.Symbols are data from Luck [69] and Fouad et al. [70] and solid lines are estimates from the association models.Luck did not publish any data for  1 ,  2 ,  3 and  4 .These values were calculated by using Luck's data for   and set of Eqs.(30).Thus, we are using Luck's published data (Pub.)and transformed data (Transf.).Colors and types of the lines refer to the same models as Figs.2-9. Figure is in a format similar to the one presented by Haghmoradi et al. [71].
values.In addition, the reported    values of Mallamace [64] show a clear maximum unlike the model results.
Finally, we calculate the fractions of k-times bonded molecules along with the fraction of free sites   at saturated conditions.For the comparisons, we use the data of Luck [69] as well as the data of Fouad et al. [70].It is worth noting that Luck did not publish any data for the fractions   (except of monomer fraction  0 ).These values were calculated by using Luck's   data and set of Eqs.(30).Fouad et al. [70] estimated the fractions of k-times bonded molecules using molecular simulations of TIP4P/2005 and iAMOEBA water models.With the exception of DAPT, all the model predictions are very close to each other.In addition, the model predictions tend to be closer to the results of TIP4P/2005 (except for DAPT).It seems that the different thermodynamic models are more consistent with the data published by Fouad et al. [70].However, one could make a case here that actually these simulation predictions by Fouad et al. [70] are inconsistent with the spectroscopy estimates of Mallamace et al. [64] and Nilsson [63].If we focus at  4 in Fig. 10, we will notice that Fouad et al. [70] estimates that  4 ≥ 0.35 at 373 K (and 1 bar).However, from Fig. 7 we see that Nilsson and Mallamace estimate that  4 ≈ 0 at 373 K (and 1 bar).

Discussion on the anomalies and the results
It is known from various studies [2,31,72] that most SAFT approaches cannot capture water's anomalous behavior.Some have been able to reproduce the density maximum by introducing a temperature dependent segment diameter to PC-SAFT [73].However, as it was demonstrated by Nezbeda et al. [31], the predicted density maximum by that version of PC-SAFT appears at all pressures, which is not correct.In addition, Marshall [26] has commented that this approach lacks theoretical rigor, especially when considering aqueous mixtures, since the theory will always force a density maximum in water.Marshall, in the development of ADPC-SAFT [24], has taken a different approach, where the segment diameter is changed depending on the amount of fully-bonded molecules ( 4 ).Marshall [24] argues that water's anomalies are caused by transitions to tetrahedral symmetry and this is one way to approach this.Even though ADPC-SAFT shows very good performance for many pure water properties, it does not predict a density maximum.
As previously stated, the origin of water's anomalous properties is still a mystery.Multiple scenarios have been proposed and explored in various studies [22].In fact, Stokely et al. [74] have shown that a microscopic cell model of water, by taking the cooperativity among Hbonds into account, is able to produce phase behaviors consistent with any of the proposed scenarios for water's phase diagram.The scenario that holds each time depends on both the amount of cooperativity as well as the strength of the directional component of the hydrogen bond.
In their studies, they combined mean field and Monte Carlo calculations which are accurate at describing water's anomalous behavior [75,76].All the scenarios in literature are connected to water's structure, phase behavior and limits of stability [22] and not explicitly to the size of the molecules.This is likely one of the reasons why the changing segment diameter methods are not that successful with water's anomalies.
Of all the scenarios the most influential one is the theory that there are two liquid structures in water and that there is a liquidliquid critical point at supercooled conditions [1,22].This theory has been supported by various computational methods [27,28,77] and even some experimental studies [23].Interestingly enough, in DAPT two structures are included in the reference fluid (hard-sphere and tetrahedral symmetry), which is, as stated earlier, quite consistent with this theory.This is probably why DAPT is more successful at predicting water's anomalous behavior compared to other models.The idea that structural fluctuations have to be included to properly capture water's behavior has also been proposed by Tsochantaris et al. [27].It is noteworthy that a similar change has been previously proposed by Marshall [25] in the development of APT.However, at that time the change affected only the dispersion term and not the association and this change was not enough to produce a density maximum.
Even if DAPT shows various anomalies, it still produces some unrealistic results such as the heat capacity maxima.While we are not certain what is the main cause of these unrealistic results, we can provide some input.As shown from Fig. 7, the amount of tetrahedrality is largely overestimated.DAPT relies heavily on this property since it affects the dispersion term as well as the association term.This is very likely one of the reasons why DAPT's performance is worse at high temperatures.If literature studies are accurate, then the amount of tetrahedrality should be insignificant or non-existent at these conditions.For  4 = 0, DAPT becomes near identical to the original thermodynamic perturbation theory and one can expect a similar qualitative behavior with the other SAFT-type models.Even with higher  4 values, DAPT shows exceptional accuracy for density at 1 bar and overall very good accuracy for density (0.98%).This high accuracy might be due to the parameterization since only vapor pressure and saturated liquid density were considered.
Regarding the  4 , studies indicate a sigmoid trend (Fig. 7).The main reason behind this trend is often not discussed, but it could be a result of steric hindrance and cooperativity (positive or negative cooperativity) effects [78,79].Due to steric effects, water molecules could be arranged in such ways that prevent them from forming 4 H-bonds.Due to cooperativity effects, the strength of hydrogen bonds could be altered in such a way that partially hydrogen bonded molecules are more stable than fully hydrogen bonded ones.This is a hypothesis of course and more studies are needed to confirm that these phenomena are responsible for the sigmoid trend.While there are some models available in literature [14], they have not been evaluated for the  4 data shown in Fig. 7.
Another issue to consider is the pressure dependence of  4 which was investigated by Tsochantaris et al. [27].They show that the pressure effect is insignificant and that higher pressures result in higher amount of tetrahedrality in SAFT-type models in contrast to literature estimates with molecular simulations.While they showed this effect for PC-SAFT, the same is true for other models as well.In Figures S1  and S2 (in Supplementary Material) we show similar results for CPA and DAPT when compared against the estimates of Pathak et al. [67].Another point to consider, which is connected to the cooperativity, is the equivalence of sites.In all cases, we have assumed that each site has the same probability to form a bond as any other.Due to this assumption, we can use Eqs.(30).How accurate is it to assume that however?By considering the various bonded fraction data of Fouad et al. [70] one can test this hypothesis by backcalculating the fractions   .If the equivalence of sites holds and Eqs.(30) are valid, then all the backcalculated fraction data should be the same.However, as shown in Figures S3 and S4, there are differences in the values.This suggests that the equivalence of sites does not fully hold.

Conclusions
Water has an extremely complex thermodynamic behavior.Many advanced thermodynamic models like PC-SAFT, CPA and SAFT-VR Mie are unable to capture these anomalies, even if they can yield on average satisfactory results.DAPT is quite successful at reproducing some anomalous properties like the density maximum.However, at high temperatures and/or high pressures DAPT shows some unrealistic results and in general is a less accurate model.At these conditions, water behaves more like a ''simple fluid'' and other models are more accurate compared to DAPT.From our analysis, SAFT-VR Mie is the most accurate, but it should be noted that there are other models in the literature like ADPC-SAFT (as shown by Nezbeda et al. [31]) that are more accurate even though they do not capture many water anomalies (like the density maximum).
This study also considered some structural properties and most notably the fraction of tetrahedrally coordinated molecules   =  4 .None of the models predicted this property accurately (qualitatively and quantitatively).This is very likely one of the reasons why DAPT fails at some conditions, since it is very dependent on this property.This shows that significant adjustments are needed to the association term.Maybe higher order perturbation is needed to account for steric hindrance and cooperativity or even an entire reformulation of the Wertheim term.

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. 3 .
Fig. 3. (a) Isothermal compressibility isobars of liquid water at different conditions.(b) Isobaric expansivity isobars of liquid water at different conditions.Lines are predicted results from the models and symbols are results from NIST[49,50] (Table1).

Fig. 5 .
Fig. 5. (a) Residual isochoric heat capacity isobars of liquid water at different conditions.(b) speed of sound isobars of liquid water at different conditions.Lines are predicted results from the models and symbols are data and results from NIST[49,50] (Table1).

Fig. 7 .
Fig. 7. Fraction of tetrahedrally coordinated molecules at 1 bar.Symbols are data from Nilsson [63] and Mallamace et al. [64].Dashed lines are estimates from molecular simulations [65-67] and solid lines are estimates from the thermodynamic models.Some published values are below the freezing point (273.15K) and they refer to supercooled liquid water at 1 bar.For the thermodynamic model estimates we considered  4 =   (Eq.(30)).
[64] used various different spectroscopy techniques (including Fourier Transformation Infrared Radiation and Raman) to characterize the relative population of water molecules that are fully hydrogen bonded, partially hydrogen bonded and not bonded.Nilsson[63] has used X-ray spectroscopy to estimate tetrahedral fractions in cooled and mostly supercooled water.They estimated the amount of tetrahedral molecules by using the O-O pair distribution function.It is discussed in the paper that the second peak in the pair distribution function is entirely caused by the tetrahedrally coordinated water molecules.Nilsson's estimates are the most recent experimental values to date.

Fig. 8 .
Fig. 8. Monomer fraction at 1 bar.Symbols are data from Mallamace et al. [64] and solid lines are estimates from the association models.

Table 1
Overview of experimental data and values from NIST used in this work.

Table 2
Parameter sets of PC-SAFT.Water is considered to form up to 4 hydrogen bonds (4C scheme).  and   are saturated pressure and saturated liquid density respectively. is the chain length (segment number),  is the temperature independent segment diameter,  is the segment energy parameter,   is the association energy and   is the association volume.

Table 3
Parameter set of CPA.In all cases, water is considered to form up to 4 hydrogen bonds (4C scheme).

Table 4
Parameter set of SAFT-VR Mie.Water is considered to form up to 4 hydrogen bonds (4C scheme).

Table 5
Parameter set of DAPT.Water is considered to form up to 4 hydrogen bonds (4C scheme).

Table 6
%AAD of all the thermodynamic models for all the thermodynamic properties considered in this work.