Thermodynamic Modeling of Calcium Sulfate Hydrates in the CaSO4–H2O System from 273.15 to 473.15 K with Extension to 548.15 K

Calcium sulfate is one of the most common inorganic salts with a high scaling potential. The solubility of calcium sulfate was modeled with the Pitzer equation at a temperature range from 273.15 to 473.15 K from published solubility data, which was critically evaluated. Only two Pitzer parameters, β(1) and β(2), with simple temperature dependency are required to model the solubility with excellent extrapolating capabilities up to 548.15 K. The stable temperature range for gypsum is 273.15–315.95 K, whereas above 315.95 K the stable phase is anhydrite. Hemihydrate is in the metastable phase in the whole temperature range, and the obtained metastable invariant temperature from gypsum to hemihydrate is 374.55 K. The obtained enthalpy and entropy changes at 298.15 K for the solubility reactions are in good agreement with literature values yielding solubility products of 2.40 × 10–05, 3.22 × 10–05, and 8.75 × 10–05 for gypsum, anhydrite, and hemihydrate, respectively. The obtained Pitzer model for the CaSO4–H2O system is capable of predicting the independent activity and osmotic coefficient data with experimental accuracy. The mean absolute average error of activity coefficient data at 298.15 K is less than 2.2%. Our model predicts the osmotic coefficient on the ice curve within 1.5% maximum error.


INTRODUCTION
Scaling or precipitation fouling, mainly forming a solid layer on equipment surfaces or piping networks, is a persistent problem encountered in many industrial processes, causing production losses, standstills, downtime and process efficiency decrease due to the reduction of equipment volume and material flow, increased heat transfer resistance, corrosion, and wearing out of construction materials. 1 Hence, scaling prevention and techniques for evaluating scaling tendencies are of considerable practical and theoretical importance in science and the engineering field.
Calcium sulfate is one of the most common inorganic salts with a high scaling potential in oil and gas production, water desalination, geothermal energy production, sulfur dioxide removal from flue gas by coal-fired power plant, 2,3 and hydrometallurgical processes of zinc, nickel, copper, and tungsten. 4−8 The demand for utilization of membrane technology is increasing and thus creating requirement for a better understanding of the solubility behavior of calcium sulfate. Moreover, the demand for process water circulation in hydrometallurgical processes will build up more and more complex and concentrated aqueous solutions, increasing the possibility of scaling. Thus, the need of thermodynamic understanding of a multicomponent aqueous solution is required, since laboratory analyses only, are not enough to comprehend the scaling potential and its variations with temperature and concentration.
In aqueous solutions, calcium sulfate forms stable hydrates with 0, 1/2, and 2 molecules of crystalline water, with the chemical names of anhydrite (AH: CaSO 4 ), hemihydrate (HH: CaSO 4 ·0.5H 2 O), and dihydrate, i.e., gypsum (DH: CaSO 4 · 2H 2 O). The stability regions of CaSO 4 hydrates depend on solution conditions, and they are influenced by temperature and composition of the aqueous solution. Therefore, understanding the phase equilibria of CaSO 4 as a function of temperature and other electrolytes is of great theoretical significance and practical importance, making it possible to estimate its scaling potential and facilitate the synthesis of calcium sulfate materials in industrial processes.
Applying the CALPHAD methodology, 9 the thermodynamic description of the binary CaSO 4 −H 2 O system is fundamental to clarify the thermodynamic behavior of calcium sulfate and its hydrates in aqueous solutions.
The aim of this study was to compile and reassess critically the experimental data of calcium sulfate and model the thermodynamic behavior of the CaSO 4 −H 2 O system up to 473.15 K. The assessment procedure was similar as used earlier for FeSO 4 −H 2 O, 10 MnSO 4 −H 2 O, 11 and NiSO 4 −H 2 O 12 systems. All experimental data used in the modeling were taken from the literature and reviewed critically. The resulting thermodynamic model was obtained using the thermodynamic equilibrium calculation program MTDATA, which uses Gibbs energy minimization routine and includes the Pitzer activity coefficient model for aqueous solutions. The CALPHAD method was used in modeling to ensure internal consistency of the thermodynamic data. 9 Furthermore, the modeling results were compared with the experimental data and other similar models to validate the accuracy of the present model and critical analysis in detail.

THERMODYNAMIC DATA
A large number of solubility measurements have been carried out for gypsum, hemihydrate, and anhydrite since the middle of 19th century. Most of the solubility data are in agreement with each other, in spite of slight deviations. However, a large controversy exists on transition temperatures between various calcium sulfate hydrates in the CaSO 4 −H 2 O system. The main reason for this problem is kinetic hindrance during phase change. Anhydrite does not crystallize with a measurable rate from water below 343.15 K, even in the presence of anhydrite seeds, and does not hydrate in several months to gypsum without gypsum seeds present. 13 Freyer and Voigt 13 reviewed the solubility of gypsum, hemihydrate, and anhydrite in the temperature range (273.15− 473.15) K at saturation pressure and pointed out that at low temperatures the stable phase is gypsum, at high temperatures it is anhydrite, whereas hemihydrate remains metastable at all temperatures. The borderlines of the solubility data points yield a transition temperature between gypsum−anhydrite from about 298.15 to 325.15 K. Within the scatter of the solubility data, the possible transition temperature of gypsum− hemihydrate phase change covers a range from less than 353.15 K to nearly 383.15 K.
On the basis of the review of Freyer and Voigt, 13

Journal of Chemical & Engineering Data
Article are 318.76, 373.85, and 471.55 K, respectively, according to the polynomial equation curves. He used a six-order polynomial for the solubility of anhydrite and four-order polynomials for hemihydrate and gypsum. He also gave several reasons for the relatively large scatter of the solubility data of calcium sulfate determined by different authors.
In this work, all experimental data were taken from the available literature. The data with less than three experimental data points in one paper or series were not considered in the modeling to obtain reasonable and reliable fitting results depending on temperature. Especially, the experimental data for anhydrite and hemihydrate under metastable conditions were critically analyzed, such as the data from 273.15 to 383.15 K for noncrystalline hemihydrate by D'Ans et al. 16 as well as data for more soluble polyformic forms such as β-anhydrite in the temperature range (278.15−383.15) K by Sborgi and Bianchi 17 and also β-anhydrite and β-hemihydrate in the temperature range (308.15−383.15) K by Power and Fabuss. 18 All of these metastable data were excluded. The solubility of gypsum measured by Farrah et al. 19 is regularly lower than other data, so it was also excluded completely from the assessment but retained for the result comparison.
All solubility data were converted to molality, mol/kg-H 2 O. The values of 136.14 and 18.015 g/mol were adopted for the molar mass of CaSO 4 and H 2 O, respectively, in the data conversion to obtain accurate data. The considered solubility data of calcium sulfate in water is collected in  15) K, respectively. All data above 473.15 K were reserved for testing the extrapolation capabilities. The criteria used to exclude or include a data point in the assessment are discussed in Section 4.1 in detail.

Pitzer Interaction
Model. The Pitzer model, one of the most widely used activity coefficient models, has been extensively used for modeling thermodynamic properties of aqueous electrolyte systems. Pitzer and his cooperators gave details of the model in the literature. 41−43 It was developed by combining the expression of Debye−Huckel electrostatic theory for long-range interactions and composition for shortrange ion-specific interactions with a virial-type expansion. Harvie and Weare 44 and Harvie et al. 45 further included unsymmetrical electrostatic mixing terms in the modified Pitzer model to improve the fit in multicomponent systems.
The mathematical expression and internal parameters of the model are clarified in eqs 1−6. Equation 1 reproduces the Debye−Huckel type contribution of the dilute solution domain. Equation 2 describes the stoichiometric mean activity coefficient (γ ± ), eq 3 formulates the osmotic coefficient (ϕ), eqs 4 and 5 give the concentration dependence of the electrolyte specific terms B ϕ and B, and eq 6 further explains the function g(x) for eq 5.
where A ϕ is the Debye−Huckel parameter, ν is the sum of the stoichiometric coefficients of cation (ν M ) and anion (ν X ), z is charge, b is an electrolyte-independent constant (b = 1.2), and the parameter values α 1 and α 2 used for 2−2 electrolyte of the Pitzer model are 1.4 and 12, respectively, in this work. Thus, the electrolyte specific parameters to be assessed are is used only for 2−2 or higher electrolytes. The C ϕ parameter is assumed to be concentration independent in the early Pitzer model version and adopted in this work. Archer 46 also suggested concentration dependency for C ϕ in a similar way to eq 4 but with different values for the internal constants. The constant values used in the above equations are the same as suggested by Pitzer 41 and also adopted by Harvie et al., 44 where concentration unit m is the molality of CaSO 4 (mol/kg of water), used throughout this paper. Activities of pure solid phases are assumed to be 1, that is, their thermodynamic properties are insensitive to pressure.

Article
The temperature dependency of Gibbs energy change of forming the solid phase according to eqs 7−9 is expressed in the following form (eq 13) The general temperature dependency of the parameters in MTDATA for the Pitzer equation (p) is 3.3. Parameter Optimization. MTDATA version 6.0 was used for parameter fitting in this work. In MTDATA, there are several excess Gibbs energy models available, including the Pitzer equation with Harvie et al. 45 modification and the NPL Pitzer model. 49 It solves thermodynamic equilibrium by using the Gibbs energy minimization technique and includes several pure substance databases and a number of excess Gibbs energy models for different kinds of solutions. It also has an assessment module to fit model parameters from experimental data. The objective function (OF) used in MTDATA is where w i is the weight of the experimental value, C i is the calculated value, E i is the experimental value, and U i is the uncertainty. All weights for the adopted experimental data, expect for rejected values, were set to 1 in the assessment. The goodness of the assessment for each experimental data point was estimated by the absolute percentage error, defined as 4. RESULTS AND DISCUSSION 4.1. Fitting Parameters. Due to the small values for solubilities of calcium sulfate hydrates in water, a new optimization approach was tested. Instead of comparing the calculated and measured molality, the difference in Gibbs energy was selected to fit the parameters of the Pitzer model. According to eq 13, at solubility limit, ΔG°(T) + RT ln K SP = 0. Thus, we obtain The uncertainty (U i ) was set to 100 J/mol for stable phases, whereas 500 J/mol was used for metastable phases. Only first three parameters A G −C G were found be adequate to describe ΔG°(T).

Journal of Chemical & Engineering Data
Article All weights for accepted experiments were set to 1 except duplicates for which a value of 0.5 was used. Several sets for temperature dependency of Pitzer parameters was tested. If all tested sets failed to model an experimental point within a given uncertainty its weight was changed to zero. However, if any of the tested sets was able to model it properly, its weight was changed to 1.
During the assessment, we discover that the parameter C ϕ has no influence on the simulation for CaSO 4 −H 2 O system as is customary for dilute solutions and was set to zero. However, it was also found out that the parameter β (0) is unnecessary, even though it was adopted in Pitzer modeling by most researchers, 21,42,44,50−53 usually with a constant value. Still, the effect of a constant value of 0.15 for β (0) used by Møller 50 and Spencer et al. 51 was tested.
The tested parameter sets and obtained objective function values are summed in Table 2, together with the assessed thermodynamic values of reactions in eqs 7−9 shown in Table  3. The four terms in Pitzer parameters are found sufficient to obtain accurate simulation results. An extra term would not improve the assessment significantly. Of the four term sets, the model C is less accurate than models D, E, J, and K according to the OF values. Set E and J have the lowest OF value of the four term sets, but their assessed ΔH°and ΔS°values for the solubility reactions are far away from the values of other sets and those obtained by HSC 9 54 (Table 3). Sets D and K have similar values for object function values, but set K produces inappropriate values of ΔH and ΔS for the gypsum solubility reaction. So the parameter set D is considered the best to model the CaSO 4 −H 2 O system and adopted in this work.
The optimized Pitzer parameters of model D obtained in this work are given in Table 4, containing the parameters collected from previous works. 21,53 The total number of fitted terms in Pitzer parameters is only four in our model compared with six terms in the model of Wang et al. 21 and nine terms in the model of Raju and Atkinson. 53 Wang et al. 21  The HMW 50,52 and SMW 51 models also used only four terms in their Pitzer parameters, but they included a neutral ion pair CaSO 4 (aq) in their model with two to three terms depending on temperature.
The relative errors, (C i − E i )/U i , between calculated (C i ) and "experimental" (E i ) values for precipitation reactions of calcium sulfate were plotted in Figure 1, in which the solid symbol means the included value whereas the hollow one means the excluded value in the assessment.
The uncertainty used for stable-phase equilibrium was 100, and 500 J/mol was used for the unstable one. Most values of the errors were close to zero. The data for gypsum and hemihydrate are in good agreement, whereas the experimental data for anhydrite deviate to some degree at high temperature. Only three anhydrite data points above temperature 373. 15

Journal of Chemical & Engineering Data
Article optimized parameters of this work. The solubility curves of each phase depending on temperature are drawn together with all collected experimental data. Predictions by recent models are also shown. The obtained difference between calculated and experimental molality for calcium sulfate hydrates is also shown in figures, where the solid symbol means the adopted value and the hollow one refers to data rejected in the assessment. The goodness of the assessment for adopted experimental data is estimated by standard deviation (SD), also known as root mean square error, defined as where i goes over all experimental points (N) with nonzero weights in the experimental set and C i is the calculated and E i is the experimental molality. The mean absolute percentage error (MAPE), also known as the absolute average relative deviation (AARD %), is used when the focus is on relative deviation The standard deviation values for fitted data are 0.00045 mol/ kg for gypsum, 0.00063 mol/kg for anhydrite, and 0.0053 mol/ kg for hemihydrate. The mean absolute percentage errors (MAPE), are 2.2, 2.5, and 7.3%, respectively. 4.2.1. Gypsum. As shown in Figure 2a,b, the calculated solubility values for gypsum are consistent with the most experimental solubility data. The temperature dependency of the solubility curve goes through most data points. Most of the differences between calculated and experimental data for gypsum are less than 0.001 mol/kg-H 2 O, with the absolute percentage error within 5%.
However, excluded points show a slight scatter, such as the data points at 273.15 K by D'Ans et al., 16 29 Still, the modeled curve around 273.15 K goes between these two rejected data points.
In addition, the excluded data from Farrah et al. 19 are scattered and lower than others. The absolute percentage errors between the calculated and experimental data for these excluded observations are above 5%.
4.2.2. Anhydrite. The calculated phase boundary for anhydrite performs excellently with the optimized parameters in the temperature range of 273.15−573.15 K, as presented in Figure 3a,b, even though the assessment of Gibbs energy according to solubility data shows high scatter at temperatures above 373.15 K in Figure 1.
The data sets with temperature ranges of 324. 65 35 were not included in the assessment. However, even though the high-temperature data was not included in the assessment, the difference between calculated and experimental data for anhydrite is less than 0.0005 mol/kg-H 2 O (Figure 3c). The extrapolating capacity over 473.15 K is also excellent (Figure 3b). The absolute percentage errors for molalities are within 5% for the adopted data.
4.2.3. Hemihydrate. Compared with gypsum and anhydrite, hemihydrate is considered as the metastable phase and attracts

Journal of Chemical & Engineering Data
Article less attention from other researchers. The assessed solubility and differences with predicted and measured solubilities are presented in Figure 4a, Table  6.
From the solubility curves calculated in this work ( Figure 5), the estimated transition temperature of gypsum to anhydrite is 315.95 K. This agrees with the research result of 315.15 ± 2 K from Azimi, 1 D'Ans, 16,57 Hill, 36 Posnjak, 58 Kelly et al., 59 Macdonald, 60 Zen, 65 Cruft and Chao, 66 Grigor'ev and Shamaev, 67 Corti and Fernandez-Prini, 68 and Kontrec et al., 69 and is in the temperature intervals determined by Freyer and Voigt 13 and Present. 14 The transition temperature of gypsum to hemihydrate has drawn less attention and it is in this study determined at 374.55 K, in agreement with Krumgalz's 15 value of 373.95 K but a little bit higher as obtained by Posnjak's. Still, the stable temperature range is 273.15−315.95 K for gypsum and above 315.95 K for anhydrite whereas hemihydrate is in a metastable phase in the whole temperature range.
The transition temperature of anhydrite to hemihydrate suggested recently by Krumgalz 15 is considered incorrect since after 471.55 K, hemihydrate would be in a stable phase instead

Journal of Chemical & Engineering Data
Article of anhydrite, which is most unlikely and has never been reported to our knowledge.   The thermodynamic properties of solubility reactions calculated from parameter values in Table 5 are collected in Table 7. As can be seen from the table, enthalpy and entropy changes calculated for the solubility reaction of gypsum by Wang et al. 21  The quality of our model is tested by comparing calculated activity coefficient data with experimental data at 298.15 K as well as activity of water on the ice curve. Neither of these data sets or similar data was used in the assessment; only solubility data was used.
The calculated activity coefficient compared with values obtained by Lilley and Briggs 73 is shown in Figure 6 with the estimated experimental error. Lilley and Briggs 73 used a value of −352.6 mV for standard electrode potential of the (Hg)Pb| PbSO 4 electrode when obtaining values for the activity coefficient. We also recalculated their results using the recent value of −352.0 ± 0.5 mV for this electrode determined by

Journal of Chemical & Engineering Data
Article Sippola and Taskinen. 74 The difference by measured and calculated activity coefficients is shown in Figure 7.
The MAPE (AARD%) values for original and recalculated data for the activity coefficient are 2.15 and 2.06%, respectively. Both Pitzer and Mayorga 43 and Rogers 55 used activity or osmotic coefficient data from the literature in their assessment of Pitzer parameters for CaSO 4 at 298.15 K. Their Pitzer parameters yield MAPE (AARD%) values for Lilley and Briggs 73 original data to 1.53 and 1.65% and for recalculated data to 1.53 and 0.93%, respectively. All of these results favor a value −352.0 mV for the standard electrode potential of the (Hg)Pb|PbSO 4 electrode.
Brown and Prue 75 measured the freezing point depression of CaSO 4 with a precision of ±0.0002 K. The measured freezing point depression and the calculated osmotic coefficient on the ice curve are put together in Table 9.
As can be seen from the table, our model predicts both measured osmotic coefficients 75 and thermodynamically estimated osmotic coefficients 76 with a standard deviation (SD) of 0.004 and a maximum error of 1.43%.
In a dilute solution, the accuracy of concentration measurements is greater. The difference of modeled water activity on the ice curve from experimental and theoretical activity of water is displayed in Figure 8. As can be seen, our model predicts the activity of water on the ice curve as better than 3 × 10 −6 .

SUMMARY AND CONCLUSIONS
The aim of this study is to give an accurate thermodynamic description of the CaSO 4 −H 2 O system and clarify its detailed thermodynamic properties for solution chemistry. The Pitzer activity coefficient approach was used to model the CaSO 4 − H 2 O system, and its parameters were assessed from critically evaluated solubility data with MTDATA software. Nine different parameter sets with varying temperature dependencies were tested. It was found that Pitzer parameter β (0) is unnecessary for modeling and its value was set to zero, even though it has been adopted in Pitzer modeling by most researchers. Thus, only Pitzer parameters β (1) and β (2) with simple temperature dependency are required to describe the CaSO 4 −H 2 O system from 273.15 to 473.15 K with good extrapolating capabilities.    The model was verified using independent activity coefficient and osmotic coefficient data not used in the assessment, which was based only on solubility data. Our model predicts the activity coefficient at 298.15 K with mean absolute percentage error (MAPE) 2.15% and activity of water on the ice curve better than 3 × 10 −6 . These results suggest that using the Gibbs energy difference as a dependent variable, the new optimizing strategy was successful.