The (apparent) Gibbs energy equation from the elements

These notes are compiled mostly following the notes from the James Connolly course at the ETH (www.perplex.ethz.ch/thermo_course), the Perple_X website (https://www.perplex.ethz.ch) and Holland and Powell (1998, Journal of Metamorphic Geology). Any error or misunderstanding is mine. Rather than a comprehensive derivation of thermodynamic quantities, these notes are intended as a cheat sheet, that could be useful for your project this week.

There is only one number (scalar) that rules all phase diagrams (assuming your variables are strictly pressure and temperature only), and I am still perplexed about the beauty resulting from it (nothing but a scalar). Although in detail the computation is rather sophisticated, the only important thing that feeds the algorithm is the Gibbs energy of the potential phases at the specific T and P that can be stable given the compositional space of your system. The trick is then to find those phases (and their proportion) having [1] the lowest Gibbs energy and [2] satisfying a mass balance conservation equation constrained by the composition of the system. So the bricks of all this construction are the Gibbs energies of the phases that can be computed by integrating the differential equation,

(1)dG=SdT+VdP

which is,

(2)G(T,P)=G0(Tr,Pr)TrTS(T,Pr)dT+PrPV(T,P)dP

where the subscript r denotes reference conditions (298.15 and 1 bar). It is implicit that the first integral is done at constant pressure (Pr) and the volume integral at a constant temperature different from Tr. Note that this expression is written differently in most common thermodynamic databases such as the ones from Holland and Powell but the expression [2] seems more natural when coming from the definition of the Gibbs function to describe isothermal, isobaric chemically closed systems (G=UTS+PV).

The computed G(T,P) values are said to be apparent because the entropy, volume and heat capacity of the elements forming the phase are not considered. The reason is that it makes computation hairy and because it is unnecessary for phase equilibria computations that need only relative differences. Remember however that the retrieved values are not absolute and therefore cannot be directly compared to other absolute computations such as those done with first principle calculations.

Reference Gibbs energy

One common source of confusion comes from different conventions about the G0(Tr,Pr) in [2] (referred to in Perple_X as SUPCRT and HSC convection). The confusion comes when even the same authors used different conventions for different versions of the datasets (such as in the case of hp02ver.dat = SUPCRT and hp11ver.dat = HSC, both by Holland and Powell). For the convection used in other datasets see https://www.perplex.ethz.ch/perplex_thermodynamic_data_file.html

SUPCRT convention

G0=Hf(Tr,Pr)Tr·S(Tr,Pr)

Remember that the enthalpy of formation assumes that the enthalpy of formation from the elements in its stable form is by convection zero (Hf=H(Tr,Pr)Helements(Tr,Pr)). The same is not true for the entropy of formation but in the SUPCRT convection, the entropy of formation of the elements is considered to be zero as well.

HSC convention

Gf(Tr,Pr)=Hf(Tr,Pr)Tr·(S(Tr,Pr)Selements(Tr,Pr))

Note that here the reference Gibbs energy is properly speaking the Gibbs energy of formation Gf because it uses Sf=S(Tr,Pr)Selements(Tr,Pr).

In most cases the experimental reference energy is obtained adiabatically and corresponds to the enthalpy of formation, this is the reason by it is tabulated as such in most datasets (usually the first column) and explains why enthalpy enters in other popular expressions of [2].

Entropy integral and heat capacity

The reversible heat transfer at constant pressure is enthalpy and can be measured experimentally at different temperatures. It is usually expressed as a polynomial function (see below).

(3)CP(T,Pr)=(HT)PrdHdT

We can relate this experimental value to entropy by noting that at constant pressure dH=TdSVdP becomes dH=TdS, and dS=dHT, therefore we can express the heat capacity in terms of entropy variations,

(4)CP=TdSdT,

thus having a way to compute S from the expression,

(5)dS=CPTdT

if a reference S(Tr,Pr) is known. Contrary to the references in thermodynamics, the reference for entropy is absolute because at 0 K the entropy of a perfect crystal is zero. This can be used to define the third law of entropy (S0(Tr,Pr)) as the integral from 0 K to Tr (298.15),

(6)0TrdS=S(Tr)S(0)=S0(Tr,Pr)=0TrCPTdT

The third law entropy can be measured by cryogenic calorimetry, but this technique is not able to measure residual configurational entropy so a more general expression for the reference entropy (and that can be derived from equilibrium experiments, as we will do this week) is

(7)S(Tr,Pr)=Sconf+0TrCPTdT

because CP is measured at ambient conditions (Pr), it is implicit that the above integration is done over a constant pressure Pr.

So the entropy integral is expressed as,

TrTS(T,Pr)dT=TrT[S(Tr,Pr)+TrTCPTdT]dT=(TTr)·STr,PrTrT[TrTCpTdT]dT

integrating by parts the double integral results in

TrTS(T,Pr)dT=T·S(Tr,Pr)+Tr·S(Tr,Pr)TTrTCpTdT+TrTCpdT

that when plugged in the main expression [2] gives

(8)G(T,P)=G0(Tr,Pr)S(Tr,Pr)T+S(Tr,Pr)TrTTrTCpTdT+TrTCpdT+PrPV(T,P)dP

noting that the reference Gibbs energy is expressed in the SUPCRT convention as

G0=Hf(Tr,Pr)Tr·S(Tr,Pr)

the term S(Tr,Pr)Tr cancel out living the innocent expression appearing in most books and early papers by Holland and Powell. In his thermo notes, J. Connolly refers to this as a perverse formulation.

(9)G(T,P)=Hf(Tr,Pr)T·S(Tr,Pr)TTrTCpTdT+TrTCpdT+PrPV(T,P)dP

I think the perversion comes from the fact that it might be not so obvious how to reach [9] from [2], specially to people like me that forgot how to make integration by parts to solve for the double entropy integral. For those like me,

(10)abuvdx=[uv]ababuvdx

so if you take u=Cp/T and v=T,

(11)TrT[TrTCpTdT]dT=TTrTCpTdTTrTCpdT

The volume term

We need still to bring the volume term to pressures higher than 1 bar. This is accomplished by using different Equations of States (EoS). Here I derive the most simple one but in practice, we will use another popular one in geosciences (Murnaghan EoS).

Isobaric expansion

It is defined as (in units of 1/K)

(12)α=1v(vT)Pr

The simplest approach is to assume that α and 1 bar is constant and independent of temperature. If so, we can solve for v (here the small letters denote specific properties, molar volume in this case). Remember also that volume needs to be expressed in energy units (J/bar/mol) so when multiplied by pressure gives Joules.

(13)αdT=dvv
(14)Tr,1T,1αdT=v(Tr,1)v(T,1)1vdv

Noting v(Tr,1)=v0, we have,

(15)α(TTr)=lnv(1,T)lnv0

Then the volume at 1 bar and T is

(16)v(T,1)=v0exp(α(TTr))

Isothermal compressibility

We can proceed similarly for the isothermal compressibility (i.e. assuming βT is a constant),

(17)βT=1v(vP)T

that is most commonly referred using its inverse, the bulk modus (KT),

(18)KT=v(Pv)T
(19)v(Tr,P)=v(Tr,1)exp(βT(PPr))

Taking both effect at the same time gives

(20)Tr,1T,1αdTTr,PrTr,PβTdP=v(Tr,Pr)v(T,P)1vdv
(21)v(T,P)=v0exp(α(TTr)βT(PPr))

Implementation in Perple_X

In Perple_X the molar Gibbs energy of a phase (g) is expressed using several equations of states. In this course, we will use one of the simplest (EoS = 2) with normal polynomials for α, CP, KT and Murnaghan type EoS for V (Holland & Powell 1998).

Heat capacity

CP(T,Pr)=c1+c2·T+c3/T2+c4·T2+c5/T+c6/T+c7/T3+c8·T3

Only c1,c2,c3,c5 are used for EoS = 2, so

(22)CP(T,Pr)=c1+c2·T+c3/T2+c5/T

Murnaghan equation for volume

Isobaric expansivity

From an experimental point of view, the isobaric expansivity is found by fitting volumetric data at 1 bar over a temperature range,

(23)α=(lnvT)Pr

That it is fitted to a polynomial of the general time

(24)α(T,Pr)=b1+b2·T+b3/T+b4/T2+b5/T

where only b1,b5,b6,b7 are used in EoS = 2,

(25)α(T,Pr)=b1+b5/T

Note that in PerpleX b5 is redundant b5=10b1,

then the volume at any temperature and 1 bar is computed as (note that we use the Taylor expansion approximation for small numbers where exp(x)1+x). This latter can be deactivated in perplex_option.dat

(26)V(T,Pr)=V0·[1+TrTα(T,Pr)dT]

Isothermal bulk modulus

The isothermal bulk modulus at the temperature T and 1bar is,

K(T,Pr)=b6+b7·(TTr)

equivalent to the HP98 expression,

K(T,Pr)=K(Tr)(11.5·104(TTr))

Note that b7 is also redundant and equal to b7=b6·1.5·104.

Volume at pressure and temperature

If the pressure derivative of the bulk modulus (K) is positive (K = 4) we retrieve the Murnaghan equation,

(27)V(T,P)=V(T,Pr)·(14P4P+K(T,Pr))1/4

by changing the notation for the bulk modulus K(T,Pr) as simply KT, and making Pr=1

(28)V(T,P)=V(T,1)·(14P4P+KT)1/4

Now, in order to compute the energy contribution due to volume at T,P we need to integrate the above expression from 1 bar to P (see Holland and Powell 1998, page 312),

(29)1PV(T,P)dP=V(T,1)KT3[(1+4PKT)3/41]

Now you have all the necessary knowledge to compute the apparent Gibbs energy of a pure mineral phase at any desired temperature and pressure conditions (EoS is expected to work at P<105 bar, <10 GPa), for higher pressure other more sophisticated EoS with thermal pressure correction effects, is needed (such as the one implemented by Holland and Powell, 2011, called Modified Tait equation of state). We will not dive into this business.