This calculator is still under development!
Please contact me regarding any issues: michael.anenburg@anu.edu.au
y axis range
y axis range
oxygen.csv
with the following content:
1000,1000
1250,1000
1500,1000
1000,10000
1250,10000
1500,10000
This interactive app allows quick and easy calculation and visualisation of various oxygen fugacity (fO2) buffers.
The calculator is fully functional. It is still work in progress and more buffers will be added in the near future.
Coming soon. In the meanwhile, a good place to start would be Frost (1991), Anenburg & O'Neill (2019) and references therein.
Coming soon, but see the two papers referenced above for more information.
Each oxygen fugacity buffer, given specified temperature and pressure conditions, constrains the chemical potential of O2 (μO2). This chemical potential is converted to fugacity, essentially expressing it as 'pseudo'-pressure, referenced to the standard state of 1 bar. It is calculated using the equation log10fO2 = (μO2° + ∫ΔVdP)/(R×T×ln10).
I am Michael Anenburg, and you can contact me at michael.anenburg@anu.edu.au.
There is no proper citation for the app as of yet, but Anenburg & O'Neill (2019) covers the overall topic of solid state fO2 buffers. Specific sources for data and thermodynamic equations used for each buffer can also be cited, as required. See 'Buffer details and references' tab for more information.
$$f\mathrm{O}_2=\left(\frac{X\mathrm{CO_2}}{X\mathrm{CO}}\right)^2\times\exp\left(\frac{\Delta_rG}{RT}\right)$$ $$\Delta_rG = 2\left(\Delta_fH_\mathrm{CO_2} + H^{\circ}_\mathrm{CO_2} - H^{\circ}_\mathrm{298.15,CO_2} - TS^{\circ}_\mathrm{CO_2}\right) - 2\left(\Delta_fH_\mathrm{CO} + H^{\circ}_\mathrm{CO} - H^{\circ}_\mathrm{298.15,CO} - TS^{\circ}_\mathrm{CO}\right) - \left(H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}\right)$$ Values for enthalpy and entropy are evaluated using a Shomate equation: $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$ All parameters are given in the NIST Chemistry WebBook “Gas phase thermochemistry data” pages for carbon dioxide, carbon monoxide, and oxygen.
This buffer is calculated at atmospheric pressures (1 bar), regardless of the pressure setting. It is only directly comparable to other buffers when they are likewise calculated at 1 bar. It is valid from 298 K to 6000 K, but does not take into account gas dissociation into monoatomic molecules.
$$f\mathrm{O}_2=\left(\frac{X\mathrm{H_2O}}{X\mathrm{H_2}}\right)^2\times\exp\left(\frac{\Delta_rG}{RT}\right)$$ $$\Delta_rG = 2\left(\Delta_fH_\mathrm{H_2O} + H^{\circ}_\mathrm{H_2O} - H^{\circ}_\mathrm{298.15,H_2O} - TS^{\circ}_\mathrm{H_2O}\right) - 2\left(\Delta_fH_\mathrm{H_2} + H^{\circ}_\mathrm{H_2} - H^{\circ}_\mathrm{298.15,H_2} - TS^{\circ}_\mathrm{H_2}\right) - \left(H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}\right)$$ Values for enthalpy and entropy are evaluated using a Shomate equation: $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$ All parameters are given in the NIST Chemistry WebBook “Gas phase thermochemistry data” pages for water, hydrogen, and oxygen.
This buffer is calculated at atmospheric pressures (1 bar), regardless of the pressure setting. It is only directly comparable to other buffers when they are likewise calculated at 1 bar. It is valid from 500 K to 6000 K, but does not take into account gas dissociation into monoatomic molecules.
$$f\mathrm{O}_2=\left(\frac{X\mathrm{SO_2}}{X\mathrm{H_2S}}\times f\mathrm{H_2O}\right)^{2/3}\times\exp\left(\frac{\Delta_rG}{RT}\right)$$ $$\Delta_rG = \frac{2}{3}\left(\Delta_fH_\mathrm{SO_2} + H^{\circ}_\mathrm{SO_2} - H^{\circ}_\mathrm{298.15,SO_2} - TS^{\circ}_\mathrm{SO_2}\right) + \frac{2}{3}\left(\Delta_fH_\mathrm{H_2O} + H^{\circ}_\mathrm{H_2O} - H^{\circ}_\mathrm{298.15,H_2O} - TS^{\circ}_\mathrm{H_2O}\right)$$ $$ - \frac{2}{3}\left(\Delta_fH_\mathrm{H_2S} + H^{\circ}_\mathrm{H_2S} - H^{\circ}_\mathrm{298.15,H_2S} - TS^{\circ}_\mathrm{H_2S}\right) - \left(H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}\right)$$ Values for enthalpy and entropy are evaluated using a Shomate equation: $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$ All parameters are given in the NIST Chemistry WebBook “Gas phase thermochemistry data” pages for water, hydrogen sulfide, and sulfur dioxide, and oxygen.
This buffer is calculated at atmospheric pressures (1 bar), regardless of the pressure setting. It is only directly comparable to other buffers when they are likewise calculated at 1 bar. It is valid from 500 K to 6000 K, but does not take into account gas dissociation into monoatomic molecules. All gases are considered as ideal gases.
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 2\left(\int_1^P V^\mathrm{NiO}\,dP - \int_1^P V^\mathrm{Ni}\,dP\right)}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill and Pownceby (1993): $$\mu\mathrm{O}_2=-478967+248.514T-9.7961T\ln{T}$$ Calibrated from 700 K to 1700 K.
Campbell et al (2009) use a Mie-Grüneisen–Birch-Murnaghan equation of state: $$P=\frac{3}{2}K_0\left(\left(\frac{V_0}{V}\right)^{7/3}-\left(\frac{V_0}{V}\right)^{5/3}\right)\left(1-\frac{3}{4}\left(4-K'_0\right)\left(\left(\frac{V_0}{V}\right)^{2/3}-1\right)\right)+\gamma_0\frac{V^{q-1}}{V_0^q}\left(E(\theta_D,T) - E(\theta_D,295)\right)$$ \(E\left(\theta_D,T\right)=3nRTD_3\left(\frac{\theta_D}{T}\right)\) where \(D_3()\) is the Debye function.
Direct calculation of \(\int_1^P V\,dP\) requires \(V=f(P)\) which is not available, so instead it is calculated using integration by parts: $$\int_{P_0}^{P_1} V\,dP = \Delta PV - \int_{V_0}^{V_1} P\,dV=P(V_1)V_1 - P(V_0)V_0 - \int_{V_0}^{V_1} P\,dV$$ $$\int P\,dV=\left(E(\theta_D,T) - E(\theta_D,295)\right)\frac{\gamma_0}{q}\left(\frac{V}{V_0}\right)^q-\frac{9K_0\left(\left(K_0'-4\right)V_0^3+\left(14-3K_0'\right)V^{2/3}V_0^{7/3}+\left(3K_0'-16\right)V^{4/3}V_0^{5/3}\right)}{16V^2}$$ Values for \(V_n\) are found numerically by solving for \(V\) which satisfy \(P(V)-P=0\).
All required parameters are given in Table 1 of Campbell et al (2009) and are calibrated to about 100 GPa and 2500 K.
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 2\left(\int_1^P V^\mathrm{FeO}\,dP - \int_1^P V^\mathrm{Fe}\,dP\right)}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill and Pownceby (1993): $$\mu\mathrm{O}_2=\begin{cases} -605568+1366.42T-182.7955T\ln{T}+0.10359T^2, & 833 < T < 1042 \\[2ex] -519113+59.129T+8.9276T\ln{T}, & 1042 \leq T \leq 1184 \\[2ex] -550915+269.106T-16.9484T\ln{T}, & 1184 < T < 1644\end{cases}$$ Calibrated at 1 bar from 833 K to 1644 K.
For FeO, Campbell et al (2009) use a Mie-Grüneisen–Birch-Murnaghan equation of state: $$P=\frac{3}{2}K_0\left(\left(\frac{V_0}{V}\right)^{7/3}-\left(\frac{V_0}{V}\right)^{5/3}\right)\left(1-\frac{3}{4}\left(4-K'_0\right)\left(\left(\frac{V_0}{V}\right)^{2/3}-1\right)\right)+\gamma_0\frac{V^{q-1}}{V_0^q}\left(E(\theta_D,T) - E(\theta_D,295)\right)$$ \(E\left(\theta_D,T\right)=3nRTD_3\left(\frac{\theta_D}{T}\right)\) where \(D_3()\) is the Debye function and \(n=2\).
Direct calculation of \(\int_1^P V\,dP\) requires \(V=f(P)\) which is not available, so instead it is calculated using integration by parts: $$\int_{P_0}^{P_1} V\,dP = \Delta PV - \int_{V_0}^{V_1} P\,dV=P(V_1)V_1 - P(V_0)V_0 - \int_{V_0}^{V_1} P\,dV$$ $$\int P\,dV=\left(E(\theta_D,T) - E(\theta_D,295)\right)\frac{\gamma_0}{q}\left(\frac{V}{V_0}\right)^q-\frac{9K_0\left(\left(K_0'-4\right)V_0^3+\left(14-3K_0'\right)V^{2/3}V_0^{7/3}+\left(3K_0'-16\right)V^{4/3}V_0^{5/3}\right)}{16V^2}$$ Values for \(V_n\) are found numerically by solving for \(V\) which satisfy \(P(V)-P=0\).
All required parameters are given in Table 1 of Campbell et al (2009) and are calibrated to about 100 GPa and 2500 K.
For Fe, the equation of state described by Dorogokupets et al (2017) is used. The equation is given in terms of the Helmholz free energy \(F\). Expressions for \(F(V,T)\) and \(P(V)\) are given in the Supplementary Files of Dorogokupets et al (2017). \(V(P,T)\) is calculated by numerically solving \(P(V)-P=0\) for the pressure of interest. \(G(P,T)\) is calculated using: $$G = F+PV$$ This Gibbs energy is used to find phase transitions, and the resulting pressures are used as the boundaries for integration. The \(\int_1^P V^\mathrm{Fe}\,dP\) integral is evaluated numerically.
This equation of state is calibrated to >300 GPa and >3000 K. The BCC α-Fe and γ-Fe phases are stable at low pressures. At high pressures, the FCC γ-Fe exists at high temperatures whereas the HCP ε-Fe is stable at lower temperatures. Evaluation of \(G\) leads to prediction of FCC γ-Fe at extreme pressures (>100–150 GPa). This phase transition is spurious and results from extrapolation beyond the calibrated range, and as such is not considered here.
Note relationship between three reactions in the Fe-O system. IW is only stable at temperature of about 830 K and higher. The invariant intersection point with the other buffers (IM and WM) shifts to lower temperature and higher oxygen fugacity with increasing pressure. See schematic diagram:
Wüstite is non-stoichiometric, and its non-stoichiometry is a function of oxygen fugacity, and is determined by the coexistence of iron as a buffering assemblage. The two studies from which wüstite data is used here contained iron, ensuring that wüstite had the appropriate stoichiometry for the IW buffer.
$$\log_{10}f\mathrm{O}_2=\log_{10}f_{\mathrm{1\,bar}}\mathrm{O}_2+\log_{10}f_P\mathrm{O}_2$$
Parameterisation at 1 bar from Fegley (2013) Table 10-16, from 950 K to 1870 K: $$\log_{10}f_{1 bar}\mathrm{O}_2=-25.7139 - \frac{19375}{T} + 11.4191\log_{10}T$$
$$\log_{10}f_P\mathrm{O}_2 = \frac{6\int_1^P V^\mathrm{hm}\,dP - 4\int_1^P V^\mathrm{mt}\,dP}{RT\ln{10}}$$ Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 7, & \mathrm{magnetite} \\ 5, & \mathrm{hematite} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
The valid pressure range of this buffer is limited by the breakdown of magnetite at high pressures: “…magnetite breaks down via a disproportionation reaction to Fe4O5 and hematite (Fe2O3)…this reaction occurs at ~9.5–11 GPa and 973–1673 K.” (Woodland et al, 2012)
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + \left(3\left(\int_1^P V^\mathrm{qz}\,dP + G^\mathrm{ex}_{P,T}\right) + 2\int_1^P V^\mathrm{mt}\,dP - 3\int_1^P V^\mathrm{fa}\,dP\right)}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill (1987):$$\mu\mathrm{O}_2=-587474+1584.427T-203.3164T\ln{T}+0.09271T^2$$ Calibrated at 1 bar from 900 K to 1420 K.
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 7, & \mathrm{fayalite\,and\,magnetite} \\ 3, & \mathrm{quartz} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
α-quartz to β-quartz transition considered using Landau theory. The transition at 1 bar is included in the expression for \(\mu\mathrm{O}_2\), hence only the pressure dependent Landau term is sought by subtracting the pressure independent Landau contribution from the total using: $$G^\mathrm{ex}_{P,T} = G^\mathrm{ex}_{P,T}(V_\mathrm{max}) - G^\mathrm{ex}_{P,T}(V_\mathrm{max}=0)$$ $$G^\mathrm{ex}_{P,T} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
Note limited stability range for the actual phase assemblage. Although the fO2 curve can be extrapolated beyond the stability range and be used as a reference curve (as ΔFMQ), it is best to avoid doing so.
Including high pressure phase transitions for Fe-ringwoodite, ferrite-type fayalite, spinel-type ferrosilite, akimotoite, and perovskite-type ferrosilite.
$$\log_{10}f\mathrm{O}_2=\frac{3\Delta G^\mathrm{Fe_2Si_2O_6} + 2\Delta G^\mathrm{Fe_3O_4} - 6\Delta G^\mathrm{Fe_2SiO_4} - \Delta G^\mathrm{O_2}}{RT\ln{10}}$$
$$\Delta G^\mathrm{solids} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{O_2} = H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}$$
\(\Delta G\) for solids is evaluated using: $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$ \(\Delta G^\mathrm{O_2}\) is evaluated using a Shomate equation and data from the NIST Chemistry WebBook “Gas phase thermochemistry data” pages. $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 7, & \mathrm{Fe_2SiO_4} \\ 5, & \mathrm{FeSiO_3} \\ 7, & \mathrm{Fe_3O_4}\end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Landau theory is used to calculate excess Gibbs energy for magnetite phase transitions: $$G^\mathrm{excess} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
Phase transitions are found by calculating \(\Delta G\) for all phases per given composition, and choosing the minimum for each combination of pressure and temperature.
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
The valid pressure range of this buffer is limited by the breakdown of magnetite at high pressures: “…magnetite breaks down via a disproportionation reaction to Fe4O5 and hematite (Fe2O3)…this reaction occurs at ~9.5–11 GPa and 973–1673 K.” (Woodland et al, 2012)
Also note complex phase relations in real high pressure assemblages containing FeO and SiO2, making FFM more of a theoretical exercise in metastable buffer curves.
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + \left(\int_1^P V^\mathrm{ReO_2}\,dP - \int_1^P V^\mathrm{Re}\,dP\right)}{RT\ln10}$$
Parameterisation at 1 bar from Pownceby and O'Neill (1994): $$\mu\mathrm{O}_2=-451020+297.595T-14.6585T\ln{T}$$ Calibrated between 850 K and 1250 K.
Pressure effect is estimated using ambient conditions molar volumes as \(\int_1^P \Delta V \,dP = 0.9917P\), and is unlikely to be reliable at pressures higher than several kbar.
Parameterisation at 1 bar from O'Neill and Nell (1997): $$\mu\mathrm{O}_2=-324563+344.151T-22.1155T\ln{T}$$ Calibrated 700 K to 1800 K.
Armstrong et al (2020) measured the properties of ruthernium metal and oxides at high temperature and pressure. For Ru metal, they use a modified Tait equation of state: $$V=V_{0,T}\left(1-a\left(1-\left(1+bP\right)^{-c}\right)\right)$$ $$\int_{P_0}^{P_1} V\,dP=V_{0,T}\left(1-a\right)\left(P_1 - P_0\right)+V_{0,T}\frac{a}{b(1-c)}\left(\left(1+bP_1\right)^{1-c} - \left(1+bP_0\right)^{1-c}\right)$$ $$V_{0,T}=V_0\exp\left(\int_{298.15}^T \alpha_0 + \alpha_1 dT\right) = V_0\exp\left(\alpha_0\left(T-298.15\right) + 0.5\alpha_1\left(T^2 - 298.15^2\right)\right)$$ $$a=\frac{1+K'_T}{1+K'_T+K_TK''_T}$$ $$b=\frac{K'_T}{K_T}-\frac{K''_T}{1+K'_T}$$ $$c=\frac{1+K'_T+K_T K''_T}{{K'_T}^2 + K'_T - K_T K''_T}$$ $$K_T=K_0+\left(\frac{\partial K}{\partial T}\right)_P\left(T-T_0\right) \quad K'_T = 4 \quad K'_T = - \frac{K'_T}{K_T}$$
Ruthenium oxide is described using a modified Tait equation of state including a thermal pressure term: $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ Assuming \(P_0=0\approx 1\,\mathrm{bar}\), then: $$\int_0^P V\,dP=V_0P\left(1-a + \frac{a}{b(1-c)P}\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1-bP\mathrm{th}\right)^{1-c}\right)\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$K'_0 = 4 \quad K''_0 = - \frac{K'_0}{K_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}\quad n=3$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Ruthenium oxide undergoes two phase transitions. It is tetragonal at low pressures, orthorhombic at intermediate pressures, and cubic at high pressures. The tetragonal to orthorhombic phase change is modelled using Landau Theory: $$G^\mathrm{ex}_{P,T} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{2}\left(Q^4 - Q^4_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = 1-\frac{T_0}{T_{c0}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \frac{T_{cP}-T}{T_{cP}} & \mathrm{otherwise} \end{cases}$$ The \(G^\mathrm{ex}_{P,T}\) term is added to the Landau-independent Gibbs energy of the tetragonal/orthogonal phase: \(G^\mathrm{total}_{P,T} = G_{P,T} + G^\mathrm{ex}_{P,T}\), where: $$G_{P,T} = \Delta_fH + \int_{T_0}^{T} Cp dT - T\left(S + \int_{T_0}^{T} \frac{Cp}{T} dT\right) $$ $$Cp = a+bT+cT^{-2}+dT^{-0.5}$$ Once the phase transition at a given \(T\) is found by solving for \(P\) which satisfies \(G_{P,T}^\mathrm{tet-orth} = G_{P,T}^\mathrm{cubic}\), oxygen fugacity is calculated by:. $$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + \left(\int_1^P V^\mathrm{RuO_2}\,dP - \int_1^P V^\mathrm{Ru}\,dP\right)}{RT\ln10}$$ If the maximum pressure is below the transition pressure, then: $$\int_1^P V^\mathrm{RuO_2}\,dP = \int_1^P V^\mathrm{RuO_2(tet-orth)}\,dP + G^\mathrm{ex}_{P,T}$$ If the maximum pressure is above the transition pressure, then: $$\int_1^P V^\mathrm{RuO_2}\,dP = \int_1^{P_\mathrm{trans}} V^\mathrm{RuO_2(tet-orth)}\,dP + \int_{P_\mathrm{trans}}^P V^\mathrm{RuO_2(cubic)}\,dP + G^\mathrm{ex}_{P_{trans},T} = $$ $$\int_1^{P_\mathrm{trans}} V^\mathrm{RuO_2(tet-orth)}\,dP + \int_1^P V^\mathrm{RuO_2(cubic)}\,dP - \int_1^{P_\mathrm{trans}} V^\mathrm{RuO_2(cubic)}\,dP + G^\mathrm{ex}_{P_{trans},T}$$
All required parameters are given in Tables 5 and 6 of Armstrong (2020). Calibrated to 25 GPa between 773 and 2500 K.
$$\log_{10}f\mathrm{O}_2=\frac{2\Delta G^\circ_\mathrm{CO} - 2\Delta G^\circ_\mathrm{C} - \Delta G^\circ_\mathrm{O_2} }{RT\ln{10}}$$
For gases, $$\Delta_rG = 2\left(\Delta_fH_\mathrm{CO} + H^{\circ}_\mathrm{CO} - H^{\circ}_\mathrm{298.15,CO} - TS^{\circ}_\mathrm{CO}\right) - \left(H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}\right)$$ Values for enthalpy and entropy are evaluated using a Shomate equation: $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$ All parameters are given in the NIST Chemistry WebBook “Gas phase thermochemistry data” pages for carbon monoxide and oxygen.
Graphite is calculated at its standard state \(\Delta G^\circ_\mathrm{C}\). Equations are described in Holland and Powell (2011) and parameters are available in their dataset version 6.33 (tc-633.txt).
$$\Delta G^\circ_\mathrm{C} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right)$$ $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$
This buffer is calculated at atmospheric pressures (1 bar), regardless of the pressure setting. It is only directly comparable to other buffers when they are likewise calculated at 1 bar. It is valid from 298 K to 6000 K, but does not take into account gas dissociation into monoatomic molecules.
$$\log_{10}f\mathrm{O}_2=\frac{\Delta G^\circ_\mathrm{CO_2} - \left(\Delta G^\circ_\mathrm{C} + \int_1^P V_\mathrm{C}\,dP\right) - \Delta G^\circ_\mathrm{O_2} }{RT\ln{10}} + \frac{\ln f\mathrm{CO}_2}{\ln{10}}$$
Oxygen fugacity values for the CCO buffer are strongly dependent on the fugacity of carbon dioxide (\(f\mathrm{CO}_2\)). There is lack of consensus between experimental studies and molecular dynamics simulations. Simulations results in consistently higher \(f\mathrm{CO}_2\) than is observed in experimental studies, which is then propagated to higher \(f\mathrm{O}_2\).
Experimental work from Frost and Wood (1997) which also summarises previous work shows that the CCO buffer follows the same slope of the FMQ buffer in (\(f\mathrm{O}_2-P\)) space, with (\(f\mathrm{O}_2\)) of the CCO buffer being roughly equal to that of FMQ at low temperatures (~800 °C), and becoming more reduced relative to FMQ at higher temperatures (ΔFMQ=–2 at ~1400 °C and 1 GPa).
In contrast, molecular dynamics simulations show a steeper slope, resulting in higher \(f\mathrm{O}_2\) values, especially at lower temperatures. For example, the simulated data of Li et al (2021) is 1 \(\log_{10}f\mathrm{O}_2\) higher than that obtained using the experimental data of Frost and Wood (1997) at 800 °C and 3 GPa. Other studies (e.g., Fu et al (2017) and Duan and Zhang (2006)) have a similar slope, although the discrepancy is of a slightly lower magnitude. The different models converge at higher temperatures, with all experimental and simulation studies agreeing to within 0.5 \(\log_{10}f\mathrm{O}_2\) unit at at 1400 °C and 3 GPa, but then deviating again at higher pressures.
As the reasons for the discrepancies are speculative without any additional studies, there is little justification to pick one model over the other. Here, two models are considered: the experimental study of Frost and Wood (1997) and the moluecular dynamics simulation of Li et al (2021). The true value is most likely lies somewhere between the two curves.
Gas species in the standard state are evaluated using a Shomate equation: $$\Delta G^\circ_\mathrm{CO_2,\,O_2} = H^{\circ} - H^{\circ}_\mathrm{298.15} - TS^{\circ}$$ $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$ All parameters are given in the NIST Chemistry WebBook “Gas phase thermochemistry data” pages for carbon dioxide, and oxygen.
Elemental carbon requires two terms: standard state \(\Delta G^\circ_\mathrm{C}\) and high pressure \(\int_1^P V_\mathrm{C}\,dP\). Equations are described in Holland and Powell (2011) and parameters are available in their dataset version 6.33 (tc-633.txt).
$$\Delta G^\circ_\mathrm{C} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + \int_1^P V\,dP$$ $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$
\(\int_1^P V_\mathrm{C}\,dP\) is calculated using a modified Tait equation of state including a thermal pressure term as described by : $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}\quad n=1$$
\(\Delta G^\circ_\mathrm{C} + \int_1^P V_\mathrm{C}\,dP\) is evaluated twice for the two carbon phases (graphite and diamond). The lower value for a given \(P\) and \(T\) is used in the expression for \(\log_{10}f\mathrm{O}_2\) given above.
Frost and Wood (1997) parameterise oxygen fugacity, with the equation provided in an erratum: $$\log_{10}f\mathrm{O}_2=4.3927-\frac{21234-117.9P+0.891P^2}{T}+7.679\times10^{-3}P+3.627\times10^{-4}P^2$$ This parameterisation is strongly overestimating \(f\mathrm{O}_2\) at low pressure (< 0.5 GPa), and underestimating \(f\mathrm{O}_2\) at high pressures (> 10 GPa). Instead, we use the parameters fitted by Frost and Wood (1997) an expanded Redlich-Kwong equation of state: $$V=\frac{RT}{P}+b-\frac{aR\sqrt{T}}{(RT+bP)(RT+2bP)}+c\sqrt{P}+dP$$ This equation is integrated to obtain: $$\ln f\mathrm{CO}_2 = \ln{P} + \frac{bP}{RT} + \frac{a\left(\ln(RT+bP)-\ln(RT+2bP)\right)}{RTb\sqrt{T}} + \frac{2cP\sqrt{P}}{3RT} + \frac{dP^2}{2RT} $$ This equation results in the correct \(f\mathrm{O}_2\) values for low pressures, and a slope consisent with simulated data at high pressures, up to about 20 GPa.
Calculated \(f\mathrm{CO}_2\) using equations of state and parameters derived from molecular dynamics simulations tend to be higher than those obtained experimentally. Here, the study of Li et al (2021) is used, as it results in the highest \(f\mathrm{O}_2\) values and can therefore be used as an upper bound for the range of simulated \(f\mathrm{CO}_2\). They used a modified Lee and Kesler equation of state: $$P = \frac{RT}{V}\left( 1+\frac{B}{V_r}+\frac{C}{V_r^2}+\frac{D}{V_r^4} + \frac{E}{V_r^5} + \frac{F}{V_r^2}\left(\beta+\frac{\gamma}{V_r^2}\right)\exp\left(-\frac{\gamma}{V_r^2}\right)\right)$$
Direct calculation of \(\int_1^P V\,dP\) requires \(V=f(P)\) which is not available, so instead it is calculated using integration by parts: $$\ln f\mathrm{CO}_2 = \int_{P_0}^{P_1} V\,dP = \Delta PV - \int_{V_0}^{V_1} P\,dV=P(V_1)V_1 - P(V_0)V_0 - \int_{V_0}^{V_1} P\,dV$$ \(\int_{V_0}^{V_1} P\,dV\) is integrated numerically, and values for \(V_n\) are found numerically by solving for \(V\) which satisfy \(P(V)-P=0\).
Including high pressure phase transitions for wadsleyite, ringwoodite, high pressure clinoenstatite, akimotoite, and bridgmanite.
$$\log_{10}f\mathrm{O}_2=\frac{\Delta G^\mathrm{MgCO_3} + \Delta G^\mathrm{MgSiO_3} - \Delta G^\mathrm{Mg_2SiO_4} - \Delta G^\mathrm{C} - \Delta G^\mathrm{O_2}}{RT\ln{10}}$$
$$\Delta G^\mathrm{solids} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{O_2} = H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}$$
\(\Delta G\) for solids is evaluated using: $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$ \(\Delta G^\mathrm{O_2}\) is evaluated using a Shomate equation and data from the NIST Chemistry WebBook “Gas phase thermochemistry data” pages. $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 7, & \mathrm{Mg_2SiO_4} \\ 5, & \mathrm{MgSiO_3} \\ 10, & \mathrm{Mg_2Si_2O_6} \\ 5, & \mathrm{MgCO_3} \\ 1, & \mathrm{C}\end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Phase transitions are found by calculating \(\Delta G\) for all phases per given composition, and choosing the minimum for each combination of pressure and temperature.
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
Including high pressure phase transitions for wadsleyite, ringwoodite, high pressure clinoenstatite, akimotoite, and bridgmanite.
$$\log_{10}f\mathrm{O}_2=\frac{\Delta G^\mathrm{CaMg(CO_3)_2} + 4\Delta G^\mathrm{MgSiO_3} - \Delta G^\mathrm{CaMgSi_2O_6} - 1\Delta G^\mathrm{Mg2SiO_4} - 2\Delta G^\mathrm{C} - 2\Delta G^\mathrm{O_2}}{RT\ln{10}}$$
$$\Delta G^\mathrm{solids} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{O_2} = H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}$$
\(\Delta G\) for solids is evaluated using: $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$ \(\Delta G^\mathrm{O_2}\) is evaluated using a Shomate equation and data from the NIST Chemistry WebBook “Gas phase thermochemistry data” pages. $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 5, & \mathrm{MgSiO_3} \\ 10, & \mathrm{Mg_2Si_2O_6} \\ 7, & \mathrm{Mg_2SiO_4} \\ 10, & \mathrm{CaMg(CO3)2} \\ 10, & \mathrm{CaMgSi2O6} \\ 1, & \mathrm{C}\end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Phase transitions are found by calculating \(\Delta G\) for all phases per given composition, and choosing the minimum for each combination of pressure and temperature.
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
Dolomite requires an additional excess energy term \(\Delta G^\mathrm{ex,dol}\) described by the symmetric formalism model. $$\Delta G^\mathrm{ex,dol} = \Delta H^\mathrm{sf} + P \Delta V^\mathrm{sf} + Q\left(W - \Delta H^\mathrm{sf} + P(W_V - \Delta V^\mathrm{sf})\right) - Q^2(W + PW_V) + RTS_\mathrm{od}$$ where \(Q\) is calculated by numerically solving $$0 = \Delta H^\mathrm{sf} + P \Delta V^\mathrm{sf} + (W + PW_V)(2Q-1) + \mathrm{fac}\frac{n}{n+1}RT\ln{\frac{n(1-Q)^2}{(1+nQ)(n+Q)}}$$ and \(S_\mathrm{od}\) is given by $$S_\mathrm{od} = \mathrm{fac}\left((1+nQ)\ln{\frac{1+nQ}{n+1}}+n(1-Q)\ln{\frac{1-Q}{n+1}}+n(1-Q)\ln{\frac{n(1-Q)}{n+1}}+n(n+Q)\ln{\frac{n+Q}{n+1}}\right)\bigg/(n+1)$$
Symmetric formalism values for dolomite are available in the Holland and Powell (2011) dataset version 6.33 (tc-633.txt): $$\Delta H^\mathrm{sf} = 11.91, \Delta V^\mathrm{sf} = 0.016, W = 11.9, W_V = 0.016, n = 1, \mathrm{fac} = 1$$
$$\log_{10}f\mathrm{O}_2=\frac{\Delta G^\mathrm{CaMg(CO_3)_2} + \Delta G^\mathrm{ex,dol} + 2\Delta G^\mathrm{SiO_2} - \Delta G^\mathrm{CaMgSi_2O_6} - 2\Delta G^\mathrm{C} - 2\Delta G^\mathrm{O_2}}{RT\ln{10}}$$
$$\Delta G^\mathrm{solids} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{O_2} = H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}$$
\(\Delta G\) for solids is evaluated using: $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$ \(\Delta G^\mathrm{O_2}\) is evaluated using a Shomate equation and data from the NIST Chemistry WebBook “Gas phase thermochemistry data” pages. $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 3, & \mathrm{SiO_2} \\ 10, & \mathrm{CaMg(CO_3)_2} \\ 10, & \mathrm{CaMgSi_2O_6} \\ 1, & \mathrm{C}\end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Phase transitions are found by calculating \(\Delta G\) for all phases per given composition, and choosing the minimum for each combination of pressure and temperature.
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
Dolomite requires an additional excess energy term \(\Delta G^\mathrm{ex,dol}\) described by the symmetric formalism model. $$\Delta G^\mathrm{ex,dol} = \Delta H^\mathrm{sf} + P \Delta V^\mathrm{sf} + Q\left(W - \Delta H^\mathrm{sf} + P(W_V - \Delta V^\mathrm{sf})\right) - Q^2(W + PW_V) + RTS_\mathrm{od}$$ where \(Q\) is calculated by numerically solving $$0 = \Delta H^\mathrm{sf} + P \Delta V^\mathrm{sf} + (W + PW_V)(2Q-1) + \mathrm{fac}\frac{n}{n+1}RT\ln{\frac{n(1-Q)^2}{(1+nQ)(n+Q)}}$$ and \(S_\mathrm{od}\) is given by $$S_\mathrm{od} = \mathrm{fac}\left((1+nQ)\ln{\frac{1+nQ}{n+1}}+n(1-Q)\ln{\frac{1-Q}{n+1}}+n(1-Q)\ln{\frac{n(1-Q)}{n+1}}+n(n+Q)\ln{\frac{n+Q}{n+1}}\right)\bigg/(n+1)$$
Symmetric formalism values for dolomite are available in the Holland and Powell (2011) dataset version 6.33 (tc-633.txt): $$\Delta H^\mathrm{sf} = 11.91, \Delta V^\mathrm{sf} = 0.016, W = 11.9, W_V = 0.016, n = 1, \mathrm{fac} = 1$$
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + \left(\int_1^P V^\mathrm{WO_2}\,dP - \int_1^P V^\mathrm{W}\,dP\right)}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill and Pownceby (1993): $$\mu\mathrm{O}_2=-569087+300.479T-15.9697T\ln{T}$$ Calibrated from 700 K to 1700 K.
Shofner et al (2016) use a Mie-Grüneisen–Birch-Murnaghan equation of state: $$P=\frac{3}{2}K_0\left(\left(\frac{V_0}{V}\right)^{7/3}-\left(\frac{V_0}{V}\right)^{5/3}\right)\left(1-\frac{3}{4}\left(4-K'_0\right)\left(\left(\frac{V_0}{V}\right)^{2/3}-1\right)\right)+\gamma_0\frac{V^{q-1}}{V_0^q}\left(E(\theta_D,T) - E(\theta_D,295)\right)$$ \(E\left(\theta_D,T\right)=3nRTD_3\left(\frac{\theta_D}{T}\right)\) where \(D_3()\) is the Debye function, and \(n=\begin{cases} 1, & \mathrm{W} \\ 3, & \mathrm{WO_2} \end{cases}\).
Direct calculation of \(\int_1^P V\,dP\) requires \(V=f(P)\) which is not available, so instead it is calculated using integration by parts: $$\int_{P_0}^{P_1} V\,dP = \Delta PV - \int_{V_0}^{V_1} P\,dV=P(V_1)V_1 - P(V_0)V_0 - \int_{V_0}^{V_1} P\,dV$$ $$\int P\,dV=\left(E(\theta_D,T) - E(\theta_D,295)\right)\frac{\gamma_0}{q}\left(\frac{V}{V_0}\right)^q-\frac{9K_0\left(\left(K_0'-4\right)V_0^3+\left(14-3K_0'\right)V^{2/3}V_0^{7/3}+\left(3K_0'-16\right)V^{4/3}V_0^{5/3}\right)}{16V^2}$$ Values for \(V_n\) are found numerically by solving for \(V\) which satisfy \(P(V)-P=0\).
\(\mathrm{WO_2}\) undergoes two phase transitions at 4 GPa and 32 GPa, and are independent of temperature in the Shofner et al (2016) model. Thus in cases where the required pressure is greater than 4 GPa or 32 GPa, the \(\int V\,dP\) integral is broken into two or three integrals, respectively, e.g.: $$\int_{P_0}^{P_1} V^{\mathrm{WO_2}}\,dP = \int_{P_0}^{4\,\mathrm{GPa}} V^{\mathrm{ap}}\,dP + \int_{4\,\mathrm{GPa}}^{32\,\mathrm{GPa}} V^{\mathrm{hp}}\,dP + \int_{32\,\mathrm{GPa}}^{P_1} V^{\mathrm{hpm}}\,dP$$
All required parameters are given in Table 1 of Shofner et al (2016) and are calibrated to about 70 GPa and 2300 K.
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 2\int_1^P V^\mathrm{cuprite}\,dP - 4\int_1^P V^\mathrm{copper}\,dP}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill (1988): $$\mu\mathrm{O}_2= -347705+246.096T - 12.9053T\ln{T}$$ Calibrated at 1 bar from 750 K to 1330 K.
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 1, & \mathrm{copper} \\ 3, & \mathrm{cuprite} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 4\int_1^P V^\mathrm{tenorite}\,dP - 2\int_1^P V^\mathrm{cuprite}\,dP}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill (1988): $$\mu\mathrm{O}_2= -292245+377.012T - 23.1976T\ln{T}$$ Calibrated at 1 bar from 800 K to 1300 K.
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 3, & \mathrm{cuprite} \\ 2, & \mathrm{tenorite} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 2\int_1^P V^\mathrm{Fe_3O_4}\,dP - 6\int_1^P V^\mathrm{FeO}\,dP}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill (1988): $$\mu\mathrm{O}_2= -581927+65.618T + 38.741T\ln{T}$$ Calibrated at 1 bar from 833 K to 1270 K.
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 7, & \mathrm{magnetite} \\ 2, & \mathrm{wustite} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt), using the ferropericlase (“fper”) end-member for wüstite.
The valid pressure range of this buffer is limited by the breakdown of magnetite at high pressures: “…magnetite breaks down via a disproportionation reaction to Fe4O5 and hematite (Fe2O3)…this reaction occurs at ~9.5–11 GPa and 973–1673 K.” (Woodland et al, 2012)
Note relationship between three reactions in the Fe-O system. WM is only stable at temperature of about 830 K and lower. The invariant intersection point with the other buffers (IM and WM) shifts to lower temperature and higher oxygen fugacity with increasing pressure. See schematic diagram:
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 0.5\int_1^P V^\mathrm{Fe_3O_4}\,dP - 1.5\int_1^P V^\mathrm{Fe}\,dP}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill (1988): $$\mu\mathrm{O}_2= -607673+1060.994T - 132.3909T\ln{T}+0.06657T^2$$ Calibrated at 1 bar from 750 K to 833 K.
For Fe, the equation of state described by Dorogokupets et al (2017) is used. The equation is given in terms of the Helmholz free energy \(F\). Expressions for \(F(V,T)\) and \(P(V)\) are given in the Supplementary Files of Dorogokupets et al (2017). \(V(P,T)\) is calculated by numerically solving \(P(V)-P=0\) for the pressure of interest. \(G(P,T)\) is calculated using: $$G = F+PV$$ This Gibbs energy is used to find phase transitions, and the resulting pressures are used as the boundaries for integration. The \(\int_1^P V^\mathrm{Fe}\,dP\) integral is evaluated numerically.
This equation of state is calibrated to >300 GPa and >3000 K. The BCC α-Fe and γ-Fe phases are stable at low pressures. At high pressures, the FCC γ-Fe exists at high temperatures whereas the HCP ε-Fe is stable at lower temperatures. Evaluation of \(G\) leads to prediction of FCC γ-Fe at extreme pressures (>100–150 GPa). This phase transition is spurious and results from extrapolation beyond the calibrated range, and as such is not considered here.
For magnetite, fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}\quad n=7$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
The valid pressure range of this buffer is limited by the breakdown of magnetite at high pressures: “…magnetite breaks down via a disproportionation reaction to Fe4O5 and hematite (Fe2O3)…this reaction occurs at ~9.5–11 GPa and 973–1673 K.” (Woodland et al, 2012)
Note relationship between three reactions in the Fe-O system. IM is only stable at temperature of about 830 K and lower. The invariant intersection point with the other buffers (IM and WM) shifts to lower temperature and higher oxygen fugacity with increasing pressure. See schematic diagram:
Calibration from Taylor and Foley (1989): $$\log_{10}f\mathrm{O}_2=11.82-\frac{33780}{T}+0.066\frac{P}{T}$$ Temperature in kelvin and pressure in bar.
Valid between 1000 °C and 1250 °C, and 9 kbar 35 kbar.
Including high pressure phase transitions for coesite, stishovite, and Fe-ringwoodite.
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + \int_1^P V^\mathrm{fayalite}\,dP - \left(\int_1^P V^\mathrm{quartz}\,dP + G^\mathrm{ex}_{P,T}\right) - 2\int_1^P V^\mathrm{Fe}\,dP}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill (1987):$$\mu\mathrm{O}_2=\begin{cases} -542941-33.182T+22.446T\ln{T}, & 900 < T < 1042 \\[2ex] -562377+103.384T+5.4771T\ln{T}, & 1042 \leq T \leq 1184 \\[2ex] -602739+369.704T-27.3443T\ln{T}, & 1184 < T < 1420\end{cases}$$ Calibrated at 1 bar from 900 K to 1420 K.
\(\int_1^P V\,dP\) for quartz and fayalite is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 7, & \mathrm{fayalite} \\ 3, & \mathrm{quartz} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
α-quartz to β-quartz transition considered using Landau theory. The transition at 1 bar is included in the expression for \(\mu\mathrm{O}_2\), hence only the pressure dependent Landau term is sought by subtracting the pressure independent Landau contribution from the total using: $$G^\mathrm{ex}_{P,T} = G^\mathrm{ex}_{P,T}(V_\mathrm{max}) - G^\mathrm{ex}_{P,T}(V_\mathrm{max}=0)$$ $$G^\mathrm{ex}_{P,T} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
All required values for quartz and fayalite can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
\(\int_1^P V\,dP\) for Fe is calculated using the equation of state described by Dorogokupets et al (2017) is used. The equation is given in terms of the Helmholz free energy \(F\). Expressions for \(F(V,T)\) and \(P(V)\) are given in the Supplementary Files of Dorogokupets et al (2017). \(V(P,T)\) is calculated by numerically solving \(P(V)-P=0\) for the pressure of interest. \(G(P,T)\) is calculated using: $$G = F+PV$$ This Gibbs energy is used to find phase transitions, and the resulting pressures are used as the boundaries for integration. The \(\int_1^P V^\mathrm{Fe}\,dP\) integral is evaluated numerically.
This equation of state is calibrated to >300 GPa and >3000 K. The BCC α-Fe and γ-Fe phases are stable at low pressures. At high pressures, the FCC γ-Fe exists at high temperatures whereas the HCP ε-Fe is stable at lower temperatures. Evaluation of \(G\) leads to prediction of FCC γ-Fe at extreme pressures (>100–150 GPa). This phase transition is spurious and results from extrapolation beyond the calibrated range, and as such is not considered here.
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 2\left(\int_1^P V^\mathrm{ilmenite}\,dP + G^\mathrm{ex}_{P,T}\right) - 2\int_1^P V^\mathrm{rutile}\,dP - 2\int_1^P V^\mathrm{Fe}\,dP}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill et al (1988):$$\mu\mathrm{O}_2=\begin{cases} -570745+97.278T+5.0278T\ln{T}, & 800 < T < 1042 \\[2ex] -578438+139.612T, & 1042 \leq T \leq 1184 \\[2ex] -574852+136.606T, & 1184 < T < 1420\end{cases}$$ Calibrated at 1 bar from 800 K to 1340 K. Note typographical error in the original paper–missing decimal point in the first expression.
\(\int_1^P V\,dP\) for rutile and ilmenite is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 3, & \mathrm{rutile} \\ 5, & \mathrm{ilmenite} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Ilmenite phase transitions considered using Landau theory. The transition at 1 bar is included in the expression for \(\mu\mathrm{O}_2\), hence only the pressure dependent Landau term is sought by subtracting the pressure independent Landau contribution from the total using: $$G^\mathrm{ex}_{P,T} = G^\mathrm{ex}_{P,T}(V_\mathrm{max}) - G^\mathrm{ex}_{P,T}(V_\mathrm{max}=0)$$ $$G^\mathrm{ex}_{P,T} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
All required values for rutile and ilmenite can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
\(\int_1^P V\,dP\) for Fe is calculated using the equation of state described by Dorogokupets et al (2017) is used. The equation is given in terms of the Helmholz free energy \(F\). Expressions for \(F(V,T)\) and \(P(V)\) are given in the Supplementary Files of Dorogokupets et al (2017). \(V(P,T)\) is calculated by numerically solving \(P(V)-P=0\) for the pressure of interest. \(G(P,T)\) is calculated using: $$G = F+PV$$ This Gibbs energy is used to find phase transitions, and the resulting pressures are used as the boundaries for integration. The \(\int_1^P V^\mathrm{Fe}\,dP\) integral is evaluated numerically.
This equation of state is calibrated to >300 GPa and >3000 K. The BCC α-Fe and γ-Fe phases are stable at low pressures. At high pressures, the FCC γ-Fe exists at high temperatures whereas the HCP ε-Fe is stable at lower temperatures. Evaluation of \(G\) leads to prediction of FCC γ-Fe at extreme pressures (>100–150 GPa). This phase transition is spurious and results from extrapolation beyond the calibrated range, and as such is not considered here.
$$\log_{10}f\mathrm{O}_2=\frac{\mu\mathrm{O}_2 + 2\left(\int_1^P V^\mathrm{usp}\,dP + G^\mathrm{ex,usp}_{T}\right) - 2\left(\int_1^P V^\mathrm{ilm}\,dP + G^\mathrm{ex,ilm}_{P,T}\right) - 2\int_1^P V^\mathrm{Fe}\,dP}{RT\ln10}$$
Parameterisation at 1 bar from O'Neill et al (1988):$$\mu\mathrm{O}_2=\begin{cases} -505563+4.99T+13.3151T\ln{T}, & 900 < T < 1042 \\[2ex] -521274+112.599T, & 1042 \leq T \leq 1184 \\[2ex] -587764+553.582T-54.37T\ln{T}, & 1184 < T < 1500\end{cases}$$ Calibrated at 1 bar from 800 K to 1340 K. Note typographical error in the original paper–missing decimal point in the first expression.
\(\int_1^P V\,dP\) for ulvöspinel and ilmenite is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 7, & \mathrm{ulvospinel} \\ 5, & \mathrm{ilmenite} \end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Ilmenite phase transitions considered using Landau theory. The transition at 1 bar is included in the expression for \(\mu\mathrm{O}_2\), hence only the pressure dependent Landau term is sought by subtracting the pressure independent Landau contribution from the total using: $$G^\mathrm{ex}_{P,T} = G^\mathrm{ex}_{P,T}(V_\mathrm{max}) - G^\mathrm{ex}_{P,T}(V_\mathrm{max}=0)$$ $$G^\mathrm{ex}_{P,T} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
All required values for ilmenite can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
Ulvöspinel requires* an additional excess energy term \(G^\mathrm{ex,usp}_{T}\) described by the symmetric formalism model. In this case, it is pressure independent, and can be simplified to: $$G^\mathrm{ex,usp}_{T} = \Delta H^\mathrm{sf} + Q(W - \Delta H^\mathrm{sf}) - Q^2 W + RTS_\mathrm{od}$$ where \(Q\) is calculated by numerically solving $$0 = \Delta H^\mathrm{sf} + W(2Q-1) + \mathrm{fac}\frac{n}{n+1}RT\ln{\frac{n(1-Q)^2}{(1+nQ)(n+Q)}}$$ and \(S_\mathrm{od}\) is given by $$S_\mathrm{od} = \mathrm{fac}\left((1+nQ)\ln{\frac{1+nQ}{n+1}}+n(1-Q)\ln{\frac{1-Q}{n+1}}+n(1-Q)\ln{\frac{n(1-Q)}{n+1}}+n(n+Q)\ln{\frac{n+Q}{n+1}}\right)\bigg/(n+1)$$
All required values for ulvöspinel can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt). Symmetric formalism values are the six last numbers for 'usp'.
*Note that symmetric formalism does not actually do anything for this buffer, because it is assumed that for atmosphereic pressure, the excess energy is already captured in the empirical \(\mu\mathrm{O}_2\) expression. Therefore, \(G^\mathrm{ex,usp}_{T,P=1}\) is subtracted from \(G^\mathrm{ex,usp}_{T,P}\), and given the pressure independence of the ulvöspinel symmetric formalism model, the two are always equal effectively leading to \(G^\mathrm{ex,usp}=0\). The derivation is presented here for completeness and compatibility with the Holland and Powell (2011) database.
\(\int_1^P V\,dP\) for Fe is calculated using the equation of state described by Dorogokupets et al (2017) is used. The equation is given in terms of the Helmholz free energy \(F\). Expressions for \(F(V,T)\) and \(P(V)\) are given in the Supplementary Files of Dorogokupets et al (2017). \(V(P,T)\) is calculated by numerically solving \(P(V)-P=0\) for the pressure of interest. \(G(P,T)\) is calculated using: $$G = F+PV$$ This Gibbs energy is used to find phase transitions, and the resulting pressures are used as the boundaries for integration. The \(\int_1^P V^\mathrm{Fe}\,dP\) integral is evaluated numerically.
This equation of state is calibrated to >300 GPa and >3000 K. The BCC α-Fe and γ-Fe phases are stable at low pressures. At high pressures, the FCC γ-Fe exists at high temperatures whereas the HCP ε-Fe is stable at lower temperatures. Evaluation of \(G\) leads to prediction of FCC γ-Fe at extreme pressures (>100–150 GPa). This phase transition is spurious and results from extrapolation beyond the calibrated range, and as such is not considered here.
The BAMM buffer (biotite-almandine-muscovite-magnetite) was suggested by Zen (1985) as a representative redox indicator for peraluminous granites and metamorphic rocks. With only pure phases considered, it lies roughly between NNO and FMQ in upper crustal conditions. Zen's original formulation consists of polynomial fits to \(\Delta G\) calculated at fixed intervals using thermodynamic data available at the time, and using constant volume changes. Here, the buffer is calculated using accurate equations of state for the solids employing the more recently available thermodynamic data.
$$\log_{10}f\mathrm{O}_2=\frac{3\Delta G^\mathrm{quartz} + 2\Delta G^\mathrm{magnetite} + \Delta G^\mathrm{muscovite} - \Delta G^\mathrm{almandine} - \Delta G^\mathrm{annite} - \Delta G^\mathrm{O_2}}{RT\ln{10}}$$
$$\Delta G^\mathrm{quartz} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{magnetite} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{muscovite} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + \int_1^P V\,dP$$ $$\Delta G^\mathrm{almandine} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + \int_1^P V\,dP$$ $$\Delta G^\mathrm{annite} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + \int_1^P V\,dP$$ $$\Delta G^\mathrm{O_2} = H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}$$
\(\Delta G\) for solids is evaluated using: $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$ \(\Delta G^\mathrm{O_2}\) is evaluated using a Shomate equation and data from the NIST Chemistry WebBook “Gas phase thermochemistry data” pages. $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 22, & \mathrm{annite} \\ 20, & \mathrm{almandine} \\ 3, & \mathrm{quartz} \\ 7, & \mathrm{magnetite} \\ 21, & \mathrm{muscovite}\end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Landau theory is used to calculate excess Gibbs energy for quartz and magnetite phase transitions: $$G^\mathrm{excess} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
At pressures below 10 GPa, the wüstite - hematite buffer is a hypothetical buffer, because the two minerals are not stable together and will react to form magnetite. It is, however, potentially useful as an indicator of the oxygen fugacity constrained by the metastable coexistence of simple end-member divalent and trivalent iron oxides (e.g. Fudali 1965). At higher pressures, magnetite is no longer stable and coexisting wüstite - hematite are stable. However, wüstite is non-stoichiomertric and data for the form of wüstite that exists in equilibrium with hematite does not exist. Instead, data for wüstite is taken from Campbell et al (2009) and is only an approximation. Real oxygen fugacity values are likely to be slightly higher than reported here.
$$\log_{10}f\mathrm{O}_2=\frac{2\Delta G^\mathrm{hematite} - 4\Delta G^\mathrm{wustite} - \Delta G^\mathrm{O_2}}{RT\ln{10}}$$
$$\Delta G^\mathrm{hematite} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{wustite} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + \int_1^P V\,dP$$ $$\Delta G^\mathrm{O_2} = H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}$$
\(\Delta G\) for solids is evaluated using: $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$ \(\Delta G^\mathrm{O_2}\) is evaluated using a Shomate equation and data from the NIST Chemistry WebBook “Gas phase thermochemistry data” pages. $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$
For hematite, fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T} \quad n=5$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Landau theory is used to calculate excess Gibbs energy for hematite phase transitions: $$G^\mathrm{excess} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
For wüstite, Campbell et al (2009) use a Mie-Grüneisen–Birch-Murnaghan equation of state: $$P=\frac{3}{2}K_0\left(\left(\frac{V_0}{V}\right)^{7/3}-\left(\frac{V_0}{V}\right)^{5/3}\right)\left(1-\frac{3}{4}\left(4-K'_0\right)\left(\left(\frac{V_0}{V}\right)^{2/3}-1\right)\right)+\gamma_0\frac{V^{q-1}}{V_0^q}\left(E(\theta_D,T) - E(\theta_D,295)\right)$$ \(E\left(\theta_D,T\right)=3nRTD_3\left(\frac{\theta_D}{T}\right)\) where \(D_3()\) is the Debye function and \(n=2\).
Direct calculation of \(\int_1^P V\,dP\) requires \(V=f(P)\) which is not available, so instead it is calculated using integration by parts: $$\int_{P_0}^{P_1} V\,dP = \Delta PV - \int_{V_0}^{V_1} P\,dV=P(V_1)V_1 - P(V_0)V_0 - \int_{V_0}^{V_1} P\,dV$$ $$\int P\,dV=\left(E(\theta_D,T) - E(\theta_D,295)\right)\frac{\gamma_0}{q}\left(\frac{V}{V_0}\right)^q-\frac{9K_0\left(\left(K_0'-4\right)V_0^3+\left(14-3K_0'\right)V^{2/3}V_0^{7/3}+\left(3K_0'-16\right)V^{4/3}V_0^{5/3}\right)}{16V^2}$$ Values for \(V_n\) are found numerically by solving for \(V\) which satisfy \(P(V)-P=0\).
All required parameters are given in Table 1 of Campbell et al (2009) and are calibrated to about 100 GPa and 2500 K.
The hypothetical FHQ buffer was suggested by Giggenbach (1987) as a “non-specific” redox buffer containing fayalite as a proxy for Fe2+ and hematite as a proxy for Fe3+. It is intended to approximate rock-buffering conditions in volcanic systems, and occasionally called the “rock” buffer. See Moretti and Stefánsson (2020) for more details.
$$\log_{10}f\mathrm{O}_2=\frac{2\Delta G^\mathrm{quartz} + 2\Delta G^\mathrm{hematite} - 2\Delta G^\mathrm{fayalite} - \Delta G^\mathrm{O_2}}{RT\ln{10}}$$
$$\Delta G^\mathrm{quartz} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{hematite} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + G^\mathrm{excess} + \int_1^P V\,dP$$ $$\Delta G^\mathrm{fayalite} = \Delta_fH^\circ + \int_{298.15}^T Cp\,dT - T\left(S^\circ + \int_{298.15}^T \frac{Cp}{T}\,dT\right) + \int_1^P V\,dP$$ $$\Delta G^\mathrm{O_2} = H^{\circ}_\mathrm{O_2} - H^{\circ}_\mathrm{298.15,O_2} - TS^{\circ}_\mathrm{O_2}$$
\(\Delta G\) for solids is evaluated using: $$\int Cp\,dT = aT+0.5bT^2-cT^{-1}+2dT^{0.5}$$ $$\int\frac{Cp}{T}\,dT = a\ln{T}+bT-0.5cT^{-2}-2dT^{-0.5}$$ \(\Delta G^\mathrm{O_2}\) is evaluated using a Shomate equation and data from the NIST Chemistry WebBook “Gas phase thermochemistry data” pages. $$H^{\circ} - H^{\circ}_\mathrm{298.15} = AT + \frac{BT^2}{2}+ \frac{CT^3}{3} + \frac{DT^4}{4} - \frac{E}{T} +F - H$$ $$S^{\circ} = A\ln{T} + BT + \frac{CT^2}{2} + \frac{DT^3}{3} - \frac{E}{2T^2} + G$$
Fugacity at pressure is calculated using a modified Tait equation of state including a thermal pressure term as described by Holland and Powell (2011): $$V=V_0\left(1-a\left(1-\left(1+b\left(P-P\mathrm{th}\right)\right)^{-c}\right)\right)$$ $$\int_{P_0}^P V\,dP=V_0\left(\left(P-P_0\right)\left(1-a\right) + \frac{a\left(\left(1+b\left(P-P\mathrm{th}\right)\right)^{1-c} - \left(1+b\left(P_0-P\mathrm{th}\right)\right)^{1-c}\right)}{b(1-c)}\right)$$ $$a=\frac{1+K'_0}{1+K'_0+K_TK''_0}$$ $$b=\frac{K'_0}{K_0}-\frac{K''_0}{1+K'_0}$$ $$c=\frac{1+K'_0+K_0 K''_0}{{K'_0}^2 + K'_0 - K_0 K''_0}$$ $$P\mathrm{th}=\alpha_0K_0\frac{\theta}{\xi_0}\left(\frac{1}{\mathrm{e}^u-1} - \frac{1}{\mathrm{e}^{u_0}-1}\right)$$ $$\theta = \frac{10636}{S/n} + 6.44 \quad\xi=\frac{u^2\mathrm{e}^u}{\left(\mathrm{e}^u-1\right)^2} \quad u=\frac{\theta}{T}$$ $$n=\begin{cases} 5, & \mathrm{hematite} \\ 7, & \mathrm{fayalite} \\ 3, & \mathrm{quartz}\end{cases}$$ $$\xi_o = \xi(T_0) \quad u_0=u(T_0) \quad T_0 = 298.15\,\mathrm{K}$$
Landau theory is used to calculate excess Gibbs energy for quartz and hematite phase transitions: $$G^\mathrm{excess} = S_\mathrm{max}\left(T_{c0}\left(Q^2+\frac{1}{3}\left(Q^6 - Q^6_0\right)\right) - T_{cP} Q^2 - T \left(Q^2_0 - Q^2 \right) \right)+P V_\mathrm{max} Q_0^2$$ $$T_{cP} = T_{c0}+P\frac{V_\mathrm{max}}{S_\mathrm{max}} \quad Q^2_0 = \sqrt{1-\frac{T_0}{T_{c0}}}$$ $$Q^2=\begin{cases} 0, & \mathrm{if}\,T > T_{cP} \\[2ex] \sqrt{\frac{T_{cP}-T}{T_{c0}}} & \mathrm{otherwise} \end{cases}$$
All required values can be obtained from the Holland and Powell (2011) dataset version 6.33 (tc-633.txt).
This model uses a combination of two equations of state. The first equation (BS), by Belonoshko and Saxena (1991) is only calibrated for pressure greater than 5 kbar. The second equation is a Redlich-Kwong (RK) equation with parameters from Connolly (2016), but it overestimates oxyxgen fugacity at pressures greater than about 50 kbar.
Here, the RK equation is used at pressures below 5 kbar and the BS equation is used at pressures above 40 kbar. In between, \(\log_{10}f\mathrm{O}_2\) is varied linearly between 5 and 40 kbar: $$\log_{10}f\mathrm{O}_2 = \log_{10}f\mathrm{O}_2^\mathrm{BS}\times x + \log_{10}f\mathrm{O}_2^\mathrm{RK}\times (1-x)$$ $$x = \frac{P - 5000}{40000 - 5000}$$
The relevant equations from Belonoshko and Saxena (1991) are $$f\mathrm{O}_2=f_\mathrm{5\,kbar}\times\exp\left(\frac{\int_5^P V\,dP}{RT}\right)$$ $$f_\mathrm{5\,kbar} = \exp\left(\frac{1000\left(f_1 + f_2T + f_3 \ln{T}\right)}{RT}\right)$$ $$P = \frac{a}{V} + \frac{b}{V^2} + \frac{c}{V^n}$$ $$a=\left(a_1 + a_2T/1000\right)\times10^4$$ $$b=\left(b_1 + b_2T/1000\right)\times10^6$$ $$c=\left(c_1 + c_2T/1000\right)\times10^9$$
Direct calculation of \(\int_1^P V\,dP\) requires \(V=f(P)\) which is not available, so instead it is calculated using integration by parts: $$\int_{P_0=\mathrm{5\,kbar}}^{P_1} V\,dP = \Delta PV - \int_{V_0}^{V_1} P\,dV=P(V_1)V_1 - P(V_0)V_0 - \int_{V_0}^{V_1} P\,dV$$ $$\int V\,dP = a\left(1-\ln V\right) + \frac{2b}{V} + \frac{cnV^{1-n}}{n-1}$$ The limits of integration are determined numerically by solving for \(V\) which satisfy \(P(V)-P=0\).
All parameters are given in Tables 4 and 5 of Belonoshko and Saxena (1991). Calibrated from 5 kbar to 1 Mbar, and from 400 K to 4000 K:
Oxygen fugacity is calculated using a Redlich-Kwong equation of state: $$f\mathrm{O}_2=\exp\left(\frac{\int_1^P V\,dP}{RT}\right)$$ $$P=\frac{RT}{V-b} - \frac{a}{\sqrt{T}V(V-b)}$$ $$a=\frac{1}{9\left(\sqrt[3]{2}-1\right)}\times\frac{R^2T_\mathrm{c}^{2.5}}{P_\mathrm{c}}\quad b=\frac{\sqrt[3]{2}-1}{3}\times\frac{RT_\mathrm{c}}{P_\mathrm{c}}$$
Direct calculation of \(\int_1^P V\,dP\) requires \(V=f(P)\) which is not available, so instead it is calculated using integration by parts: $$\int_{P_0=\mathrm{1\,bar}}^{P_1} V\,dP = \Delta PV - \int_{V_0}^{V_1} P\,dV=P(V_1)V_1 - P(V_0)V_0 - \int_{V_0}^{V_1} P\,dV$$ $$\int V\,dP = RT\left(\frac{V}{V-b} - \ln(V-b)\right) + \frac{a}{\sqrt{T}}\left(\frac{\ln{\frac{V}{V+b}}}{b} - \frac{1}{V+b}\right)$$ The limits of integration are determined numerically by solving for \(V\) which satisfy \(P(V)-P=0\). All parameters are given in Table 2 (row: “vWRK”) of Connolly (2016). Also works at pressures below 5 kbar.
This buffer is calculated at atmospheric pressures (1 bar), regardless of the pressure setting. It is only directly comparable to other buffers when they are likewise calculated at 1 bar.
Copyright © 2021 Michael Anenburg
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the ‘Software’), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED ‘AS IS’, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.