CaO-SiO 2 assessment using 3rd generation CALPHAD models

.


Introduction
Each year, 0.5 tons of cement is produced per person [1] worldwide, accounting for roughly 8 % of humanity's annual CO 2 emissions [1].If we intend to limit our CO 2 production as a species, it is vital that the high-temperature chemistry of cement production is optimised, as roughly half of the emissions arises from the raw material chemistry itself with the remainder arising from fuel usage.Producers are hesitant to explore changes to their production processes due to the enormous scales of production and tight safety, performance, and profit margins; however, change must come in order to reduce CO 2 emissions and to adapt to impurities in raw materials as available pure sources deplete.Thus, precise process models are desperately needed to accelerate these changes.
In ordinary Portland cement (OPC) clinker, the major cementitious constituents are alite and belite, which have the approximate formula Ca 3 SiO 5 and Ca 2 SiO 4 (or C 3 S and C 2 S in cement notation where C--CaO and S--SiO 2 ).The criticality of these phases in delivering the strength of cement implies that the CaO-SiO 2 binary system is also the most critical thermodynamic systems to assess for a process model of clinkerization.This is the first system in our work towards a full model for clinker thermodynamics which will include other major oxides such as Al 2 O 3 , Fe x O, Fe 2 O 3 and minor element oxides such as MgO.The solid compounds contained within this system include the four polymorphs of C 2 S (β, γ, α′ and α), three polymorphs of C 3 S (triclinic, monoclinic, and rhombohedral), the two polymorphs of CS (wollastonite and pseudowollastonite) and C 3 S 2 (rankinite) which is important for Solidia [2] type cements (inorganic binders hardening by reaction with CO 2 ).The large number of phases and their polymorphs means there is a substantial modelling challenge to overcome, as there are many reactions and transitions that link the properties of these phases.
The metals industry is another high-temperature, large-scale, and diverse production process which faces the same challenges as the cement industry.Thermodynamic researchers in this field have made enormous strides by coalescing their efforts around the so-called CAL-PHAD approach (CALculation of PHAse Diagrams for the modelling of thermodynamic properties in multicomponent systems).Key to this approach is a progressive increase in complexity, first focusing on unary phase models, such as those of a single metal, before modelling binary, ternary, and higher-order systems containing combinations of these elements.This is a critical concept, as the non-ideal interaction parameters are challenging to determine, and the lower-order models can often be directly carried over into the interaction models for the higher-order systems.Each stage carries out a thermodynamic assessment, which is a collation and critical review of all available experimental data for a particular thermodynamic system along with a careful estimate for the key thermodynamic properties (i.e., phase diagrams, enthalpic information, etc.).
The CALPHAD approach has been applied to the in CaO-SiO 2 system previously, with the first comprehensive assessment and modelling carried out by Haas et al. in 1981 [3] and it was used as the basis of all modern assessments for the CaO-SiO 2 system.While seminal and largely correct, the study by Haas et al. has typographical errors in the equations [4] which become obvious when comparing against the tabulated data.Commercial databases such as FACTSAGE [5] and SGTE also contain extensive information for the calcium silicates; however, recent updates to the CaO [6] and SiO 2 [7] end-members have not yet been incorporated and there has been a dramatic change to the assessed melting point of CaO.There is also a shortage of experimental data in the literature for the calcium silicates.For example, often the literature only includes the high temperature stable polymorph of C 3 S but other meta-stable low temperature polymorphs do exist and are very important for Portland cements [8].Additionally, the C 3 S 2 solid lacks thermodynamic data at higher temperatures, with data previously only available up to 300 K [9].Likewise, the high temperature polymorphs of C 2 S: α and α′, lack data at low temperature which we found was essential in fitting the parameters to the new CALPHAD models and to correctly calculate the full phase diagram.This work attempts to address all these issues by assimilating all available published data in its update to the CaO-SiO system.
Another aspect to improve upon previous assessments is in the choice of models used to capture the liquid phase.Previous assessments have used the two-sublattice ionic model [10,11], the modified quasichemical model [12,13], the associates model [14,15], and the cellular model [16].The main challenge in modelling the CaO-SiO system is in capturing the miscibility gap in the silica-rich oxide melt.Many previous assessments over-estimate the temperature at the onset of immiscibility due to limitations in the model chosen.Regarding extensions to higher order systems, Pelton and Blander [12] (modified quasi-chemical model) and Hillert et al. [10] (two-sublattice) were able to extend their binary models to the CaO-FeO-SiO 2 and CaO-SiO 2 -Al 2 O systems respectively, making their models promising choices for this system; however, the associates model is employed here as described in Barry et al. [15] as it is more readily transferable between thermodynamic software.The associates model also provides similar approximations for the melt phase when compared to the popular two-sublattice ionic model [17] but it is simpler due to the inclusion of oxide species on the sublattice (i.e., oxygen does not need to be considered as a separate species on a separate sublattice).This work will demonstrate the utility of the associate model.
One final improvement that can be made over the existing assessments of the CaO-SiO 2 system lies in the models for the solid phases.Previous attempts have all used polynomial models to fit the experimental data, such as the popular Shomate [18] polynomials.This approach is unreliable when extrapolating or interpolating far from experimental measurements.This makes connecting low and high temperature data difficult; however, this is vital for cements which are produced at high temperature but utilized/reacted closer to room temperature, as well as when capturing the broad range of melting temperatures in the system.The latest development in CALPHAD modelling [19] attempts to overcome these issues through the use of a so-called 3rd generation approach which makes use of Einstein terms for the low temperature data and one or more polynomial corrections to take into account the anharmonic contributions at high temperature (see Eq. ( 1)).This combination yields a single function which can be used to model data over the full temperature range from 0 K to the melting point with no junction points, thanks to the Einstein model terms capturing the underlying vibrational nature of the internal energy.In this work, all reviewed and novel experimental measurements for the solid phases in the CaO-SiO 2 system are modelled using the 3rd generation CALPHAD function, resulting in a more reliable description for the thermodynamic properties from 0 K using a single function.
The rest of this contribution is organized as follows: First, Sections and 3 describe the published data for the calcium silicates used in this study.Section 4 details the model used to describe all phases in the system and the fitting procedure for thermodynamic properties of the solid calcium silicates and the associates model, including our novel information criterion approach.The heat content measurements and DFT procedures for Ca 2 SiO 4 and Ca 3 Si 2 O 7 is detailed in Section 5. Finally, the error estimates from the resultant fit to the data points are assessed.This results in a much-needed update to the full cementitious CaO-SiO 2 assessment including a complete phase diagram.The model for this fundamental binary system will be used to extend to the CaO-Al 2 O 3 -SiO 2 -Fe x O y -MgO as part of the development of a thermodynamic database for Portland cement clinkers.
Table 2 Properties of the polymorphic transitions for all solids in the CaO-SiO 2 system found in the literature.

Literature review of the solid phases
This section details the structural properties of the individual solid compounds found within the CaO-SiO 2 system.The structural information is used to help distinguish polymorphs of individual compounds which are modelled in this work separately.Finally, the existing thermodynamic data for each structure that was used in the current study to model the solid phases is reviewed.

Structural properties
All solid compounds in this work are considered stochiometric; however, many of these phases have multiple polymorphs.For example, C 2 S has 5 polymorphs: β-C 2 S (which is a meta-stable, pressure stabilized compound which but appears in higher-order systems and in real cements), γ-C 2 S, α ′ L-C 2 S, α ′ H-C 2 S and α-C 2 S. All known polymorphs of the individual compounds are discussed here with their structural information summarised in Table 1 and the available data on the phase transitions including their respective temperatures and enthalpies presented in Table 2.
The structural properties of Ca 3 SiO 5 are complex [8].The available literature was reviewed by Regourd [20] and more recently by Dunstetter et al. [21].Ca 3 SiO 5 exists in 6 polymorphic modifications but only the high temperature rhombohedral structure is stable [22,23].The low temperature triclinic modifications (designated TI, TII and TIII in the literature [22,24]) and the intermediate monoclinic ones (designated MI and MII in the literature [25][26][27]) are all meta-stable but are observed in cement clinkers.All observed phase transformations are of the displacive type in which small shift of atoms lead to orientation disorder between SiO 4 tetrahedral units and distorted CaO 6 polyhedra.Therefore, the observed heat of transitions, as measured by quantitative thermal analysis, is small [28] (see Table 2).
Calcium orthosilicate (Ca 2 SiO 4 ) exists in five modifications.The space groups, lattice parameters, and atomic positions of all modifications are known with precision thanks to a careful study by Mumme et al. [29] using high temperature neutron diffraction combined with Rietveld analysis.The stable modification at room temperature, γ-Ca 2 SiO 4 sometimes also called calcio-olivine, transforms to α′ L -Ca 2 SiO 4 at 880 • C [30].This transformation is of the reconstructive type and accompanied by a considerable contraction of the cell volume [30].The kinetics of this transformation are very sluggish and in addition affected by a small solubility range of +0.4 wt% CaO on the lime rich side and +0.2 wt% SiO 2 on the silica rich side with respect to ideal 2:1 stoichiometry [30].α′ L -Ca 2 SiO 4 transforms into α′ H -Ca 2 SiO 4 with increasing temperature.This transition is second order and of the displacive type [29].Following this, a′ H -Ca 2 SiO 4 transforms into the high temperature α-Ca 2 SiO 4 modification above 1425 • C [30] through a semi-reconstructive transition [29].The fifth and last modification, β-Ca 2 SiO 4 (larnite), is meta-stable at atmospheric pressure, but can be fully stabilized by addition of small amounts of minor elements (e.g.0.4 % of B 2 O 3 [31]) or by quenching bulk sample from high temperature with a minimum particle size [30].The α′ L -β transformation is of the displacive type which is kinetically favoured over the stable α′ L -γ transition.The metastable transition temperature was determined by dilatometry [30] and derived from an experimental p-T phase diagram [32].
The third compound considered here, Ca 3 Si 2 O 7 , exists in two modifications: high temperature rankinite [33] and low temperature kilchoanite [34] (sometimes called Z-Phase in the literature).The transition is irreversible and takes place in the temperature interval 954-1090 • C [33].Wollastonite CaSiO 3 exists in two general modifications with a series of derived superstructures due to modified stacking sequences.The low temperature wollastonite compound [35] transforms into the pseudo-wollastonite modification [36][37][38] (sometimes also called cyclo-wollastonite in the literature).
In this work, the 4 polymorphs of C 2 S and 2 polymorphs of CS are treated as individual compounds and modelled.In addition, only rankinite-C 3 S 2 is considered for modelling in the present work due to a lack of kilchoainite date.Finally, in this assessment the polymorphs of C 3 S have been simplified to monoclinic, triclinic, and rhombohedral, as this captures the fundamental structural changes where sufficient data exists to distinguish them in this system, the modelling of the C 3 S structures is novel in this study.

Thermodynamic data
Heat capacity data for the Ca 2 SiO 4 β and γ polymorphs has been measured by Todd [39], King [9], and Grevel et al. [40] with these experiments being at low temperature allowing for the entropy at 298 K to be derived; however, no low temperature data exists for the high temperature polymorphs; α and α′.The relative enthalpy for the 4 polymorphs of Ca 2 SiO 4 measured by Coughlin [41] is used to aid fitting of the Cp for the high temperature polymorphs.Similarly, heats of formations at 298 K only exists for the low temperature polymorphs [42,43].
For the C 3 S solid, the heat capacity was also measured by Todd [39] and the relative enthalpy was measured by Gronow and Schweite [44].The heat of formation is often taken from the heat of decomposition as measured by Brunauer [45].To overcome the lack of data surrounding the monoclinic and rhombohedral Ca 3 SiO 5 polymorphs, each were given identical heat capacities derived from the triclinic form but are then individually corrected using their respective inversion temperatures and enthalpies as measured by Bigare et al. [8].
For rankinite-C 3 S 2 only the low temperature heat capacity was measured by King [2].A heat of formation measurement by Weeks [46] and EMF measurements by Benz and Wagner [47] can also be used to calculate the heat of formation and fix the absolute enthalpy; however, there is significant disagreement previous models around the Benz measurements [11], so these measurements were excluded.
The heat capacity of both polymorphs of CaSiO 3 were measured by Wagner [48] and the relative enthalpy for wollastonite was measured by Gronow [44], Courtial [49], and Pascal and Richet [50] for cyclowollastonite.The heat of formation for both polymorphs was also measured by Charlu et al. [51].For the end-members CaO and SiO 2, the latest assessments by Defferenes et al. [6] and Bajenova et al. [7] respectively, were used as they proved to be fully compatible with our present work.In the SiO 2 assessment by Bajenova et al., the tridymite polymorph was excluded as it was found not to be present in pure silica systems, this results in the use of 4 polymorphs of SiO 2 in the current work: α, β-quartz and α, β-cristobalite.Furthermore, small corrections were made for the solid CaO function due to normalisation errors in the original publication (this can be seen in the supplemental TDB file).All the reviewed experimental data used in the solid phase modelled is summarised alongside their respective uncertainties in Table A.2.

Literature review of the liquid and melt phases
The thermodynamic properties of the liquid end-members CaO and SiO 2 were recently assessed by Deffrennes et al. [6] and Bajenova et al. [7] respectively with a 3rd generation CALPHAD model, thus a reassessment of these is not required and instead their expressions, which incorporate a two-state model for the liquid phase, are adopted here.
The Phase Diagram for Ceramists [52] contains a complete study of the phase diagram of the CaO-SiO 2 system and lists the majority of the experimental liquidus points to be found for this system in the literature.Rankin and Wright [53] are the primary source and determined the melting points of C 2 S and CS as well as the CS and C 3 S 2 eutectic; However, corrections for the CaO side are required due to the dramatic change in the melting point of CaO in the recent reassessment [10].Modelling of the immiscible region of the phase diagram was most recently investigated by Tehwey and Hess [54] with similar results from Greig [55].The activities of SiO 2 and CaO in the liquid phase of the CaO-W.Abdul et al.
SiO 2 system have been extensively studied in the range of 1773-1873 K through a variety of different techniques [56][57][58][59].In this work, all available experimental invariant and liquidus data was used to model the phase diagram, the resultant parameters were then used to compute the activities and these values were compared against the experimental values to achieve a satisfactory re-calculation of the phase diagram.

3rd generation CALPHAD function
The 3rd generation CALPHAD equation used in this work, is made using a combination of Einstein terms with a small polynomial correction to allow modelling from absolute zero up to the melting point of the solid with a single function.Including zero kelvin in the fit is vital, as it is the primary method of connecting to the third law of thermodynamics and thus establishing the entropy which cannot be directly measured.Using only one expression for all temperatures also ensures thermodynamic consistency, as piece-wise expressions have discontinuities in their derivatives which may lead to unphysical discontinuities in their thermodynamic and derived properties.
The 3rd generation CALPHAD function has the following general form [19], ] . ( G α,solid (T) is the Gibbs free-energy of the solid compound α, H 0,α is a constant setting the absolute enthalpy of the compound at 0 K, R is the universal gas constant, and T is the temperature.There are N E,α Einstein terms with the ith term having an Einstein temperature θ E,i,α and prefactor D i,α .These terms have some constraints; for example, the sum of the Einstein prefactors, ∑ N i D i,α , should sum to the number of atoms in the chemical formula so that the Dulong-Petit law is approximated (i.e., lim T→∞ C v ≈ 3NR).The Einstein terms model lattice vibrations, thus the number of these terms should be correlated to the modes in the lattice; however, modern assessments such as that of CaO [6] tend to standardise around using triple-Einstein functions.The more realistic Debye model can be substituted for the Einstein model; however, it results in a non-analytic expression for the Gibbs free energy and therefore it is not a popular choice [60].Polynomial correction terms are also added with prefactors C j,α to consider the anharmonic effects at high temperature.Initially, only two non-zero prefactors exponents j = {2, 3} were recommended [19]; however, further modelling by Chen and Sundman [60] recommended using exponents j = {2, 5} to fit a wider spectrum of phases.When modelling CaO, it was found that pure oxides typically exhibit an increase in heat capacity just before melting [6] therefore, it is better to use a pair of j = {2, 13} terms or the triplet j = {2, 7, 8}.In this assessment, j = {2} is found sufficient to describe the data for all compounds.
Our modelling of the C-S system is unique as it is the first to use these 3rd generation CALPHAD functions for multicomponent oxides and a reassessment of all data from the experimental results.

Fitting methodology for the 3rd generation function
Faced with a broad spectrum of experimental data with unknown confidence intervals, a holistic approach is taken here to optimize all model parameters simultaneously.First, the stochiometric solid models are fit, then their parameters are fixed while the liquid model is fit.The liquid model fitting is carried out separately due to the relatively slow and complex calculations required, which include solving for multiphase equilibrium.
For the stochiometric phases, a custom code has been written which takes as its input experimental data in its original reported units to minimize transcription errors.This is then transformed into standardized units and constraints to carry out a non-linear least-squares regression.All experimental conditions and observations/measurements Fig. 1.An example diagnostic plot displaying the error (see Eq. ( 2)) of all experimental heat capacity measurements [9,39,40] considered while fitting C 2 S polymorphs.Outer horizontal lines indicate the 2σ deviation limit that should contain 95 % of the data set, while the dashed horizontal lines indicate the 1σ deviation which should contain 2/3 of data outside of 1.These guidelines are held true by adjusting the error estimates reported here.
have uncertainties attached to them, and these uncertainties are propagated through to the final uncertainty estimate for that measurement which is then used to scale its contribution to the error term.To illustrate this, consider a single experimental observation/measurement enumerated by α for a thermodynamic property (e.g., a measured enthalpy of reaction, or heat capacity), y obs,α with estimated measurement uncertainty σ y,α .This is observed at experimental conditions (i.e., temperature and pressure), x → obs,α , which may also have associated uncertainties σ → x,α .The corresponding value predicted by the model is given by y ) where c → are the fit coefficients.The reduced error, e α , is then given by the following expression, where the total uncertainty, σ α , assuming no correlations between variables and Gaussian distributed noise, can be estimated using the standard expression for uncertainty propagation, The uncertainty estimates on the observed value σ y typically account for all sources of uncertainty in the experimental conditions, σ → x ; however, transition temperatures are entered as reactions with an "observed" Gibbs free energy change of zero as this is mathematically exact at a transition.To capture the uncertainty of the measurement in terms of free energy, uncertainty is assigned to the transition temperature, which is much more natural and easily propagated to the freeenergy by this approach.
To reduce human bias, the estimation of the uncertainties of each dataset is optimised deterministically.In the first instance the author's stated uncertainties for the experimental measurement are entered.After an initial fit, the deviation between the model predictions and the experimental measurements are used to scale the uncertainty estimates.The error limits suggested by Haas et al. [3] are employed, where the uncertainties are adjusted with the goal to keep 95 % of the data within 2σ (standard deviations) of the final fit and no more than a 1/3 of data outside of 1σ, to prevent over and under weighting of datasets.This method is a deterministic way of estimating uncertainty where there is conflicting data.
The approach for determining the uncertainties has been discussed in the literature [61]; however, there is no solid consensus on the best approach as of yet.The simple approach used here still requires that data which is clearly unphysical or seriously conflicting is excluded entirely as it can result in non-convergence.But even with this limitation, the approach applied here is a powerful method that allows corroboration of results from separate authors and techniques.Fig. 1 gives an example error plot of the heat capacity data for the C 2 S polymorphs used within this assessment, error plots like this are produced for each compound in this system to compare and validate the uncertainties placed on the differing experimental results.A table of all experimental data alongside their initial and final uncertainties used in our fit to achieve acceptable results (reduced error ~1) can be seen in Table A.2.To provide an initial estimate of the coefficients for the global optimisation, individual phases are added to the optimisation sequentially as determined by the required data dependencies before a global optimization is carried out.The underlying minimiser used is the COIN-OR IPOpt solver [62] and all derivatives required for the minimisation, uncertainty propagation, and thermodynamic functions are generated using a custom Computer Algebra System (CAS).This allows arbitrary thermodynamic functions to be fit with little work and facilitated early comparisons with Debye functions, as well as transformations into other thermodynamic calculation software.

Thermodynamic model for the phase diagram
The melt model used here approximates the melt as a mixture of neutral species, called the "associates".This model is popular in systems where the melt phase exhibits short range order [63] which is the case in metal-oxygen system such as C-S and is readily implemented in popular thermodynamic software.The fictitious associates are formed from the liquid end-members but may be normalised in different ways; for example, in previous studies of the C-S system, authors have normalised associates to have 2 non-oxygen atoms per mole of species [14].In this assessment a simple direct linear mixing of the end members with corrections for the enthalpy and entropy were used.This was chosen as it produced satisfactory results and it is similar to the normalisation found in MTDATA thermodynamic database [64].The Gibbs free energy of the associates is given by the following expression, where ΔH CxSy,0 and TΔS CxSy,0 are constant enthalpy and entropy correction terms respectively for the intermediate associates and are zero for the pure end-member associates.
The set of possible associates, and thus potential values of x and y in Eq. ( 4), is determined by considering the four-way bonding of silicon to other atoms via an oxygen species.A bridging oxygen links two Si atoms, while and non-bridging oxygens link Si to a metal atom.This leads to the concept of Q n species, where n is the number of non-bridging oxygens atoms: The Q 2 and Q 0 species combined with the metal atom (Ca 2+ ) to form neutral species gives the C 2 S and CS associates which are used here.These associates appear to be justified by experimental spectroscopic evidence for their stable existence in significant quantities in the melt [65,66].In addition, the congruent melting points at these compositions is more easily captured by using associates at these compositions.
The overall expression for the free energy of the associate model for the melt phase is then given by the following expression, where x α is the mole fraction of the associate, α.The final term of Eq. ( 5) is a Redlich-Kister (RK) correction for non-ideal melt interactions (excess Gibbs energy) between the associates and L α,β,k is the coefficient for the kth-order correction term between associate α and β, where and k ≥ 0. The interaction terms may be composed of several corrections, where h α,β,k is an enthalpy shift term, s α,β,k T is an entropy shift term, and c p,α,β,k TlnT is a heat-capacity shift term.In this work, the final term was not necessary to fit melt interactions parameters.The fitting and choice of the associates and RK parameters are discussed in the following section, where a novel melt fitting approach is also presented.

Akaike information criterion approach to modelling phase diagram
When modelling the melt phase, there are an enormous number of choices that can be made.Fixing on the associates model, the choice remains on which Q-species associates to include, in addition to the arbitrarily large number of possible RK terms between the associates.The choice of which RK terms and associates to include follows some basic rules.The minimum number of terms required to gain a good fit is used, and lower-order interaction terms, i.e., small k, are preferred over higher-order terms as they are less likely to induce spurious features.Finally, the values of the coefficients are also bounded by experience.For example, the term s α,β,k directly affects the entropy of the phase at absolute zero, thus it's use can introduce a thermodynamic inconsistency and odd behaviour at low temperatures.While liquid phases in our W. Abdul et al. application are rarely extrapolated to very low temperatures this may happen numerically during the solution to equilibrium leading to spurious re-entrant liquid phases.In addition, the s α,β,k term becomes more important in the free-energy at high temperatures, thus positive values may lead to fictitious miscibility gaps appearing.The C-S system has a real miscibility gap between CS and S (see Fig. 8), so positive L α,β,k terms are required between pairs of associates spanning this region (i.e., L CS,S,k , L C2S,S,k , and/or L C,S,k ); however, care must be taken that this is the only location where they appear.Even with this guidance, assessments typically proceed by trial and error, with researchers adding terms to address deficiencies in the calculated diagram after successive calculations.While experts can produce high quality diagrams this way in short order, it is not deterministic, and it is difficult to know if the optimal choice has been made.Furthermore, for higher-order fitting to ternary systems such as CaO-Al 2 O 3 -SiO 2 or CaO-MgO-SiO 2 the problem of term selection and optimisation becomes factorially harder.
Here a novel automated mathematical optimization approach is used to select the RK and associate terms in the liquid to find the best combination of terms.To prevent overfitting, the objective function that is optimised is the corrected Akaike information criteria [67], where AiCc the corrected Akaike information criteria, k is the number of parameters (in this case the number of h α,β,k and s α,β,k terms), n is the number of experimental data points, and L is the maximum of the likelihood function when the h α,β,k and s α,β,k terms are optimised to fit the experimental data.The log likelihood is the sum of the error terms in Eq.
(2).The corrected equation is used to avoid overfitting in cases where the number of samples is sparse.Using the AICc, the model with the lowest score is the model that fits the data the best when traded off against the complexity of the model.This approach is entirely inspired by ESPEI [68] where it is currently used for thermochemical data fitting (i.e.stochiometric solids).
The method proceeds from an initial parameterisation (i.e., an ideal system) in a stage-by-stage manner.At each stage, all zero interaction terms are individually enabled and optimised with the current set of enabled terms.All RK terms up to an order of k = 3 are tested as higherorder terms are avoided where possible.At the end of each stage, the interaction term that is enabled permanently is the term that decreases the AiCc parameter the most.These stages continue until enabling any additional parameter no longer decreases the AICc.The optimisation of all interaction parameters against the experimental data is carried out using ESPEI.ESPEI uses a Bayesian approach that finds the maximum likelihood at each stage (more details on this technique can be found in [68]).Unfortunately, this approach currently struggles to optimize miscibility gaps when the initial guess does not capture the experimental data points well.Thus, for this assessment the optimisation for the C-S system is begun by including h CS,S,0 , h CS,S,1 , h CS,S,2 , s CS,S,0 s CS,S,1 and s CS,S,2 with an initial hand-fit.Anecdotally, the automatic approach appears to function well in systems where there is no miscibility gap, and the issue appears to only be one of local/global convergence, thus we expect it to work with further improvements to ESPEI.From this initial point, the algorithm outlined here adds a h C,C2S,1 then adds a s C2S,S,0 term.At this point, any additional terms do not lower the AiCc score (the full parameter set can be seen in Table 8).

Experimental and computational methods
The data that is available for this system is discussed in Section 2. There are several key measurements which should be repeated, or missing data that must be ascertained to complete a thermodynamic assessment targeted at modelling cement clinker formulation.This section details the preparation and synthesis methods used to obtain heat content measurements for C 2 S and C 3 S 2 .To supplement the experimental datasets, DFT calculations were also performed to generate additional theoretical data.This is important especially for high temperature phases like α′-C 2 S, for which low temperature data does not exist (these low temperature data is important in the fitting of the 3rd generation function).

Synthesis of C 3 S 2 and γ-C 2 S
The starting materials are CaCO 3 (Alpha Aesar, 99.5 %) and SiO 2 (Sibelco, 99 %) which are dried in a muffle furnace at 500 • C in air for 24 h.Rankinite (Ca 3 Si 2 O 7 ) is prepared from a dry powder mixture of CaCO 3 (6.25 ± 0.01 g) and SiO 2 (2.50 ± 0.01 g).The mixture is homogenized in a mechanical mixer for 24 h and pelletized using a hydraulic workshop press (4.5 tons on 3 g of powder).The first thermal treatment is carried out at 1000 • C/1H to calcine the limestone (CaCO 3 ) and is followed by a high temperature isothermal annealing step at 1380 • C for 12 h.The annealing process is repeated four times at 1420 • C/12 h with intermediate grinding steps between each thermal treatment.The compositional analysis reveals a loss on ignition of 31.9444%.For Ca 2 SiO 4 , a powder mixture of CaCO 3 (6.97± 0.01 g) and SiO 2 (2.09 ± 0.01 g) is calcined at 1000 • C for 1 h and annealed at 1550 • C for 8 h followed by slow cooling to room temperature.A loss on ignition of 34.0844 % is recorded.The mineralogical composition of the samples was carried out using X-ray diffraction (X'Pert Pro MPD, PANalytical), with CuKα radiation of λ = 1.5419Å.

Heat content measurements
A MHTC96 (Setaram) high temperature calorimeter with a drop measuring head (Type S, 1300 • C max) was used for the heat content measurements.The measurements are carried out using an alumina crucible with a platinum inlay (diameter 13 mm and height 46 mm).The measurements were performed under flowing Ar atmosphere.The calorimeter is calibrated using pure alumina (Al 2 O 3 , 99.95 %) pieces which are dried and calcined at 1400 • C for 12 h.The pellets have a mass ranging from 15 to 150 mg (uncertainty ±0.02 mg).Three series of 6 or 7 alumina samples are dropped for each temperature at 600 • C, 700 • C, 800 • C, 900 • C, 1000 • C, 1100 • C, and 1200 • C and the associated heat effect recorded and integrated using the Calypso software package (Setaram).The resulting calibration constant at a given temperature is calculated using heat content values for Al 2 O 3 from the FTOxid database [5].The measured heat flux area (in μV⋅s) for each sample is plotted as a function of heat content of Al 2 O 3 .The calibration constant (in μV⋅s/J) was calculated by performing a linear regression passing through origin for each series at a given temperature.The uncertainty of the calibration constant corresponds to twice the standard variation of the slope.The temperature of introduction was measured with a thermocouple (uncertainty ±0.1 K).The variation of the high temperature calorimeter is ±0.4K and a return to the baseline is observed within 25 min.
For the heat content measurement of Ca 3 Si 2 O 7 and Ca 2 SiO 4 , pellets are produced with sample weights ranging from 15 to 80 mg (uncertainty ±0.02 mg).Again, two series of 6/7 pellet are dropped at each temperature of calibration.The observed integrated heat effects (surface area of peak in μV⋅s) are plotted as a function of the phase amount introduced into the calorimeter.A linear regression forced passing through the origin is performed to determine the slope.The uncertainty of the measurement corresponds to twice the uncertainty on the slope (95 % confidence interval).The measured heat content corresponds to the slope divided by the calibration constant.Again, the introduction temperature was measured by thermocouple (±0.1 K) and corrected to 298.15 K using the low temperature heat capacity data from King [9] for Ca 3 Si 2 O 7 and Ca 2 SiO 4 .
Additionally, the melting point of Ca 3 Si 2 O 7 was measured by differential scanning calorimetry in a Labsys™ Evo apparatus designed by SETARAM.Calibration of the DSC equipment was carried out using the melting point of pure NaCl (99.9 % Alfa Aesar) as well as the α-β W. Abdul et al. transition and the melting point of K 2 SO 4 (99.5 % Alfa Aesar) as reference samples with 5 K and 10 K heating and cooling rates.

DFT calculations
The ground state properties of all experimentally stable calcium silicates were calculated using density functional theory (DFT) [69,70].The Vienna ab-initio software package (VASP v. 5.4.4) was used for the theoretical ground state calculations [71,72].The semi-local SCAN exchange functional [73] was used for the many-body exchange-correlation interaction.The GW PAW potentials (Ca_sv_gw, Si_gw and O_gw_new in VASP notation) were used for the calculations as they give the most accurate description of the lattice volume of SiO 2 in its a-quartz structure at low temperature.For Ca, the 3s, 3p and 4s orbitals, for Si, the 3s and 3p orbitals and for O the 2s and 2p orbitals were considered as valence states in the calculations.The cut-off energy was set to 800 eV.
Finite temperature properties such as the Helmholtz free energy F and heat capacities at constant volume/pressure C V /C P are approximated using lattice dynamics theory [76].The phonon spectrum of selected compounds is determined using the frozen phonon (supercell) method.The vibrational modes are calculated using the phonopy code [77] coupled to VASP.The convergence criteria for the Hellman-Feynman forces is set to 10 − 8 eV/Å to avoid residual strain in the lattice.To treat the long-range interaction of the macroscopic electrical field induced by polarization of collective motions near the Γ-point, the Born effective charges Ζ for the independent atoms in the primitive cell are calculated and used to calculate the non-analytical term corrections to the dynamical matrix.

Results and discussion
The following sections discuss the results of the experiments, fitting, and overall phase diagram in comparison to the literature.

Calorimetric measurements
The results of the heat content measurements for Rankinite are summarised in Table 3.No experimental data exists in the literature on the heat capacity or heat content of Rankinite at these high temperatures, thus this data is compared to the estimated data of Haas et al. [3] and the commercial FTPs database [5].The estimated heat content together with the experimental values from 298.15 K to 1471.65 K are presented in Fig. 2.There is some agreement with the estimation of Haas but the calculated values using FTPs clearly overestimate the heat content compared with the experimental results.
The DSC results for C 3 S 2 were only interpretable on heating due to the peritectic nature of its melting reaction.On cooling, the sample  The obtained results compared to calculated values using the FTPs database are plotted in Fig. 3.
Above 1072 K, γ-and α′-Ca 2 SiO 4 measurements do not agree with calculated values using FTPs.For α′-Ca 2 SiO 4 , our measured values are 24 to 32 kJ lower.However, there is good agreement with Haas' [3] data and Coughlin's [43] values.The measured values for γ-Ca 2 SiO 4 in the temperature range 972-1072 K are in agreement with the calculated ones from FTPs [5] but also those of Haas [3] and Coughlin [43].

DFT calculations
The ground state energies of all considered compounds are This work_Fit ( Fig. 3. γ, α′-Ca 2 SiO 4 Heat contents measurement.Lines were calculated using data from FactSage FTPs data [5] and Haas [3]: γ-Ca 2 SiO 4 300 K to 1072.0 K, α′-Ca 2 SiO 4 1171.39K to 1471.05K. Coughlin [43] values for α′-Ca 2 SiO 4 have been corrected with enthalpy of transition γ → β from Forest [34] to have the same reference state.summarised in Table 5 together with the calculated transition energies at 0 K for the various Ca 2 SiO 4 and CaSiO 3 modifications.All data are given in eV per formula unit.For Ca 3 SiO 5 , only the monoclinic modification is calculated as it is the only modification with completely occupied lattice sites.For Ca 2 SiO 4 , the γ, β and α′ modifications are calculated.The high temperature α-modification has a more complicated crystal structure with partially occupied lattice sites and could therefore not be calculated with sufficient precision.The calculated heat of formation at 0 K is presented in Table 6.The data for Ca 2 SiO 4 have an additional contribution from the Zero Point energy calculations (− 0.4 kJ mol − 1 for γ-Ca 2 SiO 4 and − 2.5 kJ mol − 1 for β-Ca 2 SiO 4 and α′-Ca 2 SiO 4 ).

Thermodynamic modelling
The resultant coefficients for the solid phases in the C-S system are given in the Table A.1.Additionally, a TDB file for this system has also been provided and has been tested for use in Pandat [78] and ESPEI [68] A table showing the full set of experimental data along with their reported and calculated uncertainties is given in Table A.2.These fits allow for an evaluation of the new data for C 3 S 2 and α ′ and α C 2 S. A comparison of the H 298 and S 298 values between the new fit and the values of Haas et al. [3] is given in Table 7.Most of the values are within the expected error of each other; however, disagreement exists for the H 298 of rankinite-Ca 3 Si 2 O 7 .This is expected and was previously reported by Hillert et al. [11].The H 298 and S 298 of α′-Ca 2 SiO 4 also disagree with the Haas values, this is explained with the use of new data for both the enthalpy and heat capacity determined in this study.Fig. 4 shows the Fig. 4. Graph of enthalpy of formations at 298 K of the stable phases at room temperature (C 3 S is metastable at room temperature), relative to the single oxides compared against the experimental data from Ref [42,43,46,79].DFT data is the calculated enthalpy of formation at 0 K in the present work.

Table 8
Table showing

Table 9
Summary of transitions and invariant reactions in the C-S system from the calculated phase diagram compared to experimental data in parentheses.The composition values correspond to the liquid phase.calculated enthalpy of formation relative to the oxides for the compounds within this system (for which experimental data exists) and compares them against the experimental data.There is generally a good agreement between these values with the only discrepancy being the C 3 S-triclinic phase (x = 0.25).This is due to the incorporation of the value of Brunauer et al. [45] who They measured the enthalpy of formation of C 3 S from the C + C 2 S reaction and provides some conflict to the value of King [42].
The resulting optimised liquid phase model including the Gibbs energy terms of the associates is presented in Table 8 and the calculated phase diagram is shown in Fig. 8.The calculated liquidus shows good agreement with the available experimental data.The miscibility gap is hard to reproduce as was noted by previous assessments [10,16], but is well produced in this assessment with only a slight disagreement on the SiO 2 rich side.To facilitate precise reading of key points of the phase diagram a summary of the invariant points is given in Table 9.While activity data was not fit in this assessment, the calculated model activities of SiO 2 and CaO alongside experimental comparisons can be seen in Figs. 5 and 6 respectively.The calculated values of SiO 2 activity is in agreement with the experimental values of Baird and Taylor [58], with excellent agreement at high CaO concentrations across all experimental values [56,58,59].Additionally, the calculated CaO activities are in line with the results of previous models of the CaO-SiO 2 system [10,13], with good agreement against the values of Carter and Macfarlane [80].
Direct comparisons of the model reported here can be made with the assessment conducted by Besmann and Spear [14] who also made use of an associates model.They used an additional C 3 S associate, which is unfavourable as it is not evidenced to exist via spectroscopic measurements [65,66].Furthermore, adding an associate adds additional complexity from the resulting corrections as well possible complications in the higher order system which is why it was avoided in this assessment.Our agreement with the experimental data demonstrates this associate is not needed from a pure fitting perspective.Not enough information is given to compare the parameterisation further; however, we estimate they use a total of 10 parameters to capture the liquid, which can be compared to our 12-parameter fit.Our assessment uses two more parameters to capture the miscibility gap, thus this may be the only difference.Hillert et al. carried out two assessments for the CaO-SiO 2 system using two variations on an ionic two-sublattice model [10,11], which due to the species they chose to represent the melt phase is a comparable model to our associates model (i.e.(Ca 2+ ), (SiO 4 ).In one model they included SiO 3 − 2 while in the other it was removed.Their models use 6 interaction terms between their sublattice species and SiO 2 to capture the liquid, whereas our approach shows not all the species require interaction parameters between them to fit the phase diagram, additionally they use higher-order ideal interactions (i.e., k = 3 RK terms) in the previously preferred model.The liquid model can also be compared against the spectroscopic data of Mysen et al. [66].Point equilibrium calculations to obtain the abundance of the associate species in our melt phase and the results are as follows are presented in Fig. 7.In the first instance, the large abundances of the C 2 S and CS species experimentally determined in the melt phase, justify the use of both associates in the final model.Furthermore, considering the large uncertainties reported by Mysen et al. [66] for these measurements, the calculated C 2 S values agree with measurement values.The calculated CS values are seen not to match the experimental values, which may be due to the large interaction parameter sets between the CS and SiO 2 associates which is required to produce the miscibility gap.Finally, Mysen et al. reported the presence of C 3 S 2 and CS 2 species but they were found in such low quantities (>10 %), they were removed from the modelling, though this difference could account for some of the disagreement with the experimental data.

Conclusion
The CaO-SiO 2 is the most important system for Portland cement clinkers but modelling this system is challenging due to the miscibility gap, the large number of phases, and the polymorphs requiring revaluation, as well as the diversity of thermodynamic data.Novel results for the relative enthalpy of rankinite Ca 3 Si 2 O 7 at high temperatures, the relative enthalpy of γ-Ca 2 SiO 4 including the transition to the α′ modifications, and phonon calculations to determine the low temperature heat capacity of α′ Ca 2 SiO 4 are presented here.This has led to a full reevaluation of the solid compounds in this system, as well as the inclusion of updates to the CaO and SiO 2 end-members [6,7] and the inclusion of the Ca 3 SiO 5 polymorphs within the model, an important feature when extending the system.The produced model appears to accurately reproduce all known experimental data, including all solid and liquid phases.Further work will focus on extending this assessment to carry out equilibrium calculations like those done by Hanein et al. [4,84] for Portland cement clinkers.The key improvement being the inclusion of precipitation of phases from the melt which is a key dynamic in that application.
A modern fitting approach has been developed and tested in this paper with the use of the 3rd generation CALPHAD functions and the associate's model.The use of information criterion to select parameters appears to have been successful, although this is only a preliminary result in the development of these techniques.Further assessments to complete the cement clinker system will provide ample opportunity for testing.These techniques will allow extensions with relative ease and work is already underway to complete the CaO-Al2O 3 -SiO 2 -Fe 2 O 3 which is a system that contains all of the major clinker phases found in Portland cement, with the iron-oxide being a key component in correctly capturing the bulk melting point.Completion of this modelling will allow us to carry prediction calculations for the composition of the clinker, which is the ultimate goal of our work.

Funding
This work was supported by Nanocem as part of the CP17 project.

Declaration of competing interest
Wahab Abdul reports financial support was provided by Nanocem.a These measurements were not fit and are only reported for comparison, their fitted uncertainty refers to the uncertainty to ensure the measurement is 1 sigma from the fit.

Fig. 7 .
Fig. 7. Abundance (Mole Frac) of C 2 S and CS species in the melt compared against spectroscopic data from Mysen et al. (including experimental uncertainties) [66].

Table 1
Structural information, polymorphism & transition data for the solid phases in the CaO-SiO 2 binary system as found in the literature.

Table 4 (
Ca 2 SiO 4 heat contents by an inversed drop calorimetric technique.The temperature of transition γ → α′ is situated near 1120 K. Therefore, in this present work the measured heat contents values include the γ-Ca 2 SiO 4 to α′-Ca 2 SiO 4 in the temperature range 972.86-1471.65K. [52]′) Ca 2 SiO 4 heat contents measurements from 972.86 to 1471.05K.Transition from γ to α′ occurs at 1170 K. contained γ-Ca 2 SiO 4 and CS which did not completely convert to C 3 S 2 .However, the melting temperature obtained is 1737 ± 5 K, which is identical to the value reported by Levin[52].The results of the heat content measurements for belite are summed up in Table 4.The heat content measurements of γ → α′ Ca 2 SiO 4 have not been available previously.Coughlin & O'Brien [41] measured the β, α′-

Table 5
Calculated ground state energies at 0 K in eV per formula unit.

Table 6
Calculated heat of formation at 0 K in kJ mol − 1 relative to the oxides CaO and α-SiO 2 (low quartz) and calculated heat of transitions for the different Ca 2 SiO 4 and CaSiO 3 modifications.

Table 7 H
[3]and S 298 of the compounds in the C-S system as calculated by the model presented here and compared to values from Haas et al.[3].
a This value is extrapolated beyond the temperature limits of the curve.W.Abdul et al.

Table A . 1
(1)le of the optimised coefficients for the solid phases in the C-S system, for use in the 3rd generation CALPHAD function given in Eq.(1).Table of all experimental measurements from literature including the uncertainty used in the fitting.The reported uncertainties are the uncertainties reported by the respective authors for that measurement and the fitted uncertainties are uncertainties used in the simultaneous fit.In the case where the uncertainties were not stated or incorrect by the author, standard error estimates of 5 % were used depending on the experiment type.Additionally, where no conflicting data exists, we have not presented a "fitted uncertainty" (see discussion in Section 3).All values are in SI units.
(continued on next page) W.Abdul et al.