Mathematical Background

Here is a bit of background on the methods used to calculate thermoelastic and thermodynamic properties in BurnMan. More detail can be found in the cited papers.

Endmember Properties

Calculating Thermoelastic Properties

To calculate the bulk (\(K\)) modulus, shear modulus (\(G\)) and density (\(\rho\)) of a material at a given pressure (\(P\)) and temperature (\(T\)), optionally defined by a geotherm) and determine the seismic velocities (\(V_S, V_P, V_\Phi\)), one uses an Equation of State (EoS). Currently the following EoSs are supported in BurnMan:

  • Birch-Murnaghan finite-strain EoS (excludes temperature effects, [Poi91]),
  • Birch-Murnaghan finite-strain EoS with a Mie-Grüneisen-Debye thermal correction, as formulated by [SLB05].
  • Birch-Murnaghan finite-strain EoS with a Mie-Grüneisen-Debye thermal correction, as formulated by [MBR+07].
  • Modified Tait EoS (excludes temperature effects, [HuangChow74]),
  • Modified Tait EoS with a pseudo-Einstein model for thermal corrections, as formulated by [HollandPowell11].
  • Compensated-Redlich-Kwong for fluids, as formulated by [HP91].

To calculate these thermoelastic parameters, the EoS requires the user to input the pressure, temperature, and the phases and their molar fractions. These inputs and outputs are further discussed in User input.

Birch-Murnaghan (isothermal)

The Birch-Murnaghan equation is an isothermal Eulerian finite-strain EoS relating pressure and volume. The negative finite-strain (or compression) is defined as

(1)\[f=\frac{1}{2}\left[\left(\frac{V}{V_0}\right)^{-2/3}-1\right],\]

where \(V\) is the volume at a given pressure and \(V_0\) is the volume at a reference state (\(P = 10^5\) Pa , \(T\) = 300 K). The pressure and elastic moduli are derived from a third-order Taylor expansion of Helmholtz free energy in \(f\) and evaluating the appropriate volume and strain derivatives (e.g., [Poi91]). For an isotropic material one obtains for the pressure, isothermal bulk modulus, and shear modulus:

(2)\[P=3K_0f\left(1+2f\right)^{5/2}\left[1+\frac{3}{2}\left(K_0^\prime -4\right) f\right],\]
(3)\[\begin{split}K_{T}=(1+2f)^{5/2}\biggl[ & K_0+(3K_0{K}^\prime_{0}-5K_0)f\\ &+\frac{27}{2}(K_0{K}^\prime_{0}-4K_0)f^2 \biggr],\end{split}\]
(4)\[\begin{split}G = (1+& 2f)^{5/2} \biggl[G_0+(3K_0{G}^\prime_{0}-5G_0)f\\ & +(6K_0{G}^\prime_{0}-24K_0-14G_{0} + \frac{9}{2}K_{0}{K}^\prime_{0})f^2 \biggr].\end{split}\]

Here \(K_0\) and \(G_0\) are the reference bulk modulus and shear modulus and \(K_0^\prime\) and \({G}^\prime_{0}\) are the derivative of the respective moduli with respect to pressure.

BurnMan has the option to use the second-order expansion for shear modulus by dropping the \(f^2\) terms in these equations (as is sometimes done for experimental fits or EoS modeling).

Modified Tait (isothermal)

The Modified Tait equation of state was developed by [HuangChow74]. It has the considerable benefit of allowing volume to be expressed as a function of pressure. It performs very well to pressures and temperatures relevant to the deep Earth [HollandPowell11].

(5)\[\begin{split}\frac{V_{P, T}}{V_{1 bar, 298 K}} &= 1 - a(1-(1+bP)^{-c}), \\ a &= \frac{1 + K_0'}{1 + K_0' + K_0K_0''}, \\ b &= \frac{K_0'}{K_0} - \frac{K_0''}{1 + K_0'}, \\ c &= \frac{1 + K_0' + K_0K_0''}{K_0'^2 + K_0' - K_0K_0''}\end{split}\]

Mie-Grüneisen-Debye (thermal correction to Birch-Murnaghan)

The Debye model for the Helmholtz free energy can be written as follows [MBR+07]

\[\begin{split}\mathcal{F} &= \frac{9nRT}{V}\frac{1}{x^3} \int_{0}^{x} \xi^2 \ln (1-e^{-\xi}) d\xi, \\ x &= \theta / T, \\ \theta &= \theta_0 \exp \left( \frac{\gamma_0-\gamma}{q_0} \right), \\ \gamma &= \gamma_0 \left( \frac{V}{V_0} \right)^{q_0}\end{split}\]

where \(\theta\) is the Debye temperature and \(\gamma\) is the Grüneisen parameter.

Using thermodynamic relations we can derive equations for the thermal pressure and bulk modulus

\[\begin{split}P_{th}(V,T) &= - \frac{\partial \mathcal{F(V, T)}}{\partial V}, \\ &= \frac{3 n \gamma R T}{V} D(x), \\ K_{th}(V,T) &= -V \frac{\partial P(V, T)}{\partial V}, \\ &= \frac{3 n \gamma R T}{V} \gamma \left[ (1-q_0 - 3 \gamma) D(x) + 3\gamma \frac{x}{e^x - 1} \right], \\ D(x) &= \frac{3}{x^3} \int_{0}^{x} \frac{\xi^3}{e^{\xi} - 1} d\xi\end{split}\]

The thermal shear correction used in BurnMan was developed by [HamaSuito98]

\[G_{th}(V,T) &= \frac{3}{5} \left[ K_{th} (V, T) - 2\frac{3nRT}{V}\gamma D(x) \right]\]

The total pressure, bulk and shear moduli can be calculated from the following sums

\[\begin{split}P(V, T) &= P_{\textrm{ref}}(V, T_0) + P_{th}(V, T) - P_{th}(V, T_0), \\ K(V, T) &= K_{\textrm{ref}}(V, T_0) + K_{th}(V, T) - K_{th}(V, T_0), \\ G(V, T) &= G_{\textrm{ref}}(V, T_0) + G_{th}(V, T) - G_{th}(V, T_0)\end{split}\]

This equation of state is substantially the same as that in SLB2005 (see below). The primary differences are in the thermal correction to the shear modulus and in the volume dependences of the Debye temperature and the Gruneisen parameter.

HP2011 (thermal correction to Modified Tait)

The thermal pressure can be incorporated into the Modified Tait equation of state, replacing \(P\) with \(P-\left(P_{\textrm{th}} - P_{\textrm{th0}}\right)\) in Equation (5) [HollandPowell11]. Thermal pressure is calculated using a Mie-Grüneisen equation of state and an Einstein model for heat capacity, even though the Einstein model is not actually used for the heat capacity when calculating the enthalpy and entropy (see following section).

\[\begin{split}P_{\textrm{th}} &= \frac{\alpha_0 K_0 E_{\textrm{th}} }{C_{V0}}, \\ E_{\textrm{th}} &= 3 n R \Theta \left(0.5 + \frac{1}{ \exp(\frac{\Theta}{T}) - 1 }\right), \\ C_{V} &= 3 n R \frac{(\frac{\Theta}{T})^2\exp(\frac{\Theta}{T})}{(\exp(\frac{\Theta}{T})-1)^2}\end{split}\]

\(\Theta\) is the Einstein temperature of the crystal in Kelvin, approximated for a substance \(i\) with \(n_i\) atoms in the unit formula and a molar entropy \(S_i\) using the empirical formula

\[\Theta_i=\frac{10636}{S_i/n_i + 6.44}\]

SLB2005 (for solids, thermal)

Thermal corrections for pressure, and isothermal bulk modulus and shear modulus are derived from the Mie-Grüneisen-Debye EoS with the quasi-harmonic approximation. Here we adopt the formalism of [SLB05] where these corrections are added to equations (2)(4):

(6)\[\begin{split}P_{th}(V,T) &={\frac{\gamma \Delta \mathcal{U}}{V}}, \\ K_{th}(V,T) &=(\gamma +1-q)\frac{\gamma \Delta \mathcal{U}}{V} -\gamma ^{2} \frac{\Delta(C_{V}T)}{V} ,\\ G_{th}(V,T) &= -\frac{\eta_{S} \Delta \mathcal{U}}{V}.\end{split}\]

The \(\Delta\) refers to the difference in the relevant quantity from the reference temperature (300 K). \(\gamma\) is the Grüneisen parameter, \(q\) is the logarithmic volume derivative of the Grüneisen parameter, \(\eta_{S}\) is the shear strain derivative of the Grüneisen parameter, \(C_V\) is the heat capacity at constant volume, and \(\mathcal{U}\) is the internal energy at temperature \(T\). \(C_V\) and \(\mathcal{U}\) are calculated using the Debye model for vibrational energy of a lattice. These quantities are calculated as follows:

\[\begin{split}C_V &= 9nR\left ( \frac{T}{\theta}\right )^3\int_{0}^{\frac{\theta}{T}}\frac{e^{\tau}\tau^{4}}{(e^{\tau}-1)^2}d\tau, \\ \mathcal{U} &= 9nRT\left ( \frac{T}{\theta} \right )^3\int_{0}^{\frac{\theta}{T}}\frac{\tau^3}{(e^{\tau}-1)}d\tau, \\ \gamma &= \frac{1}{6}\frac{\nu_{0}^2}{\nu^{2}}(2f+1)\left [ a_{ii}^{(1)} +a_{iikk}^{(2)}f\right ], \\ q &= \frac{1}{9\gamma}\left [ 18\gamma^{2}-6\gamma -\frac{1}{2} \frac{\nu^{2}_0}{\nu^2}(2f+1)^{2}a_{iikk}^{(2)} \right ], \\ \eta_S &=-\gamma-\frac{1}{2}\frac{\nu_{0}^2}{\nu^2}(2f+1)^{2}a_{S}^{(2)}, \\ \frac{\nu^2}{\nu^2_0} &= 1+a_{ii}^{(1)}f+\frac{1}{2}a_{iikk}^{(2)}f^2, \\ a_{ii}^{(1)} &= 6\gamma _0, \\ a_{iikk}^{(2)} &= -12\gamma _0+36\gamma_{0}^{2}-18q_{0}\gamma_0, \\ a_{S}^{(2)} &=-2\gamma _0-2\eta_{S0},\end{split}\]

where \(\theta\) is the Debye temperature of the mineral, \(\nu\) is the frequency of vibrational modes for the mineral, \(n\) is the number of atoms per formula unit (e.g. 2 for periclase, 5 for perovskite), and \(R\) is the gas constant. Under the approximation that the vibrational frequencies behave the same under strain, we may identify \(\nu/\nu_0 = \theta/\theta_0\). The quantities \(\gamma_0\), \(\eta_{S0}\) \(q_0\), and \(\theta_0\) are the experimentally determined values for those parameters at the reference state.

Due to the fact that a planetary mantle is rarely isothermal along a geotherm, it is more appropriate to use the adiabatic bulk modulus \(K_S\) instead of \(K_T\), which is calculated using

(7)\[K_S=K_{T}(1+\gamma \alpha T),\]

where \(\alpha\) is the coefficient of thermal expansion:

(8)\[\alpha=\frac{\gamma C_{V}V}{K_T}.\]

There is no difference between the isothermal and adiabatic shear moduli for an isotropic solid. All together this makes an eleven parameter EoS model, which is summarized in the Table below. For more details on the EoS, we refer readers to [SLB05].

User Input Symbol Definition Units
V_0 \(V_{0}\)
Volume at P = \(10^5\)
Pa , T = 300 K
m \(^{3}\) mol \(^{-1}\)
K_0 \(K_{0}\) Isothermal bulk modulus at P=10^5 Pa, T = 300 K Pa
Kprime_0 \(K^\prime_0\) Pressure derivative of \(K_{0}\)  
G_0 \(G_{0}\) Shear modulus at P = \(10^5\) Pa, T = 300 K Pa
Gprime_0 \(G^\prime_0\) Pressure derivative of \(G_{0}\)  
molar_mass \(\mu\) mass per mole formula unit kg \(\mathrm{mol}^{-1}\)
n n number of atoms per formula unit  
Debye_0 \(\theta_{0}\) Debye Temperature K
grueneisen_0 \(\gamma_{0}\) Grüneisen parameter at P = \(10^5\) Pa, T = 300 K  
q0 \(q_{0}\) Logarithmic volume derivative of the Grüneisen parameter  
eta_s_0 \(\eta_{S0}\) Shear strain derivative of the Grüneisen parameter  

This equation of state is substantially the same as that of the Mie-Gruneisen-Debye (see above). The primary differences are in the thermal correction to the shear modulus and in the volume dependences of the Debye temperature and the Gruneisen parameter.

Compensated-Redlich-Kwong (for fluids, thermal)

The CORK equation of state [HP91] is a simple virial-type extension to the modified Redlich-Kwong (MRK) equation of state. It was designed to compensate for the tendency of the MRK equation of state to overestimate volumes at high pressures and accommodate the volume behaviour of coexisting gas and liquid phases along the saturation curve.

\[\begin{split}V &= \frac{RT}{P} + c_1 - \frac{c_0 R T^{0.5}}{(RT + c_1 P)(RT + 2 c_1 P)} + c_2 P^{0.5} + c_3 P, \\ c_0 &= c_{0,0} T_c^{2.5}/P_c + c_{0,1} T_c^{1.5}/P_c T, \\ c_1 &= c_{1,0} T_c/P_c, \\ c_2 &= c_{2,0} T_c/P_c^{1.5} + c_{2,1}/P_c^{1.5} T, \\ c_3 &= c_{3,0} T_c/P_c^2 + c_{3,1}/P_c^2 T\end{split}\]

Calculating Thermodynamic Properties

So far, we have concentrated on the thermoelastic properties of minerals. There are, however, additional thermodynamic properties which are required to describe the thermal properties such as the energy, entropy and heat capacity. These properties are related by the following expressions:

(9)\[\mathcal{G}= \mathcal{E} - T\mathcal{S} + PV = \mathcal{H} - TS = \mathcal{F} + PV\]

where \(P\) is the pressure, \(T\) is the temperature and \(\mathcal{E}\), \(\mathcal{F}\), \(\mathcal{H}\), \(\mathcal{S}\) and \(V\) are the molar internal energy, Helmholtz free energy, enthalpy, entropy and volume respectively.

HP2011

(10)\[\begin{split}\mathcal{G}(P,T) &= \mathcal{H}_{\textrm{1 bar, T}} - T\mathcal{S}_{\textrm{1 bar, T}} + \int_{\textrm{1 bar}}^P V(P,T)~dP, \\ \mathcal{H}_{\textrm{1 bar, T}} &= \Delta_f\mathcal{H}_{\textrm{1 bar, 298 K}} + \int_{298}^T C_P~dT, \\ \mathcal{S}_{\textrm{1 bar, T}} &= \mathcal{S}_{\textrm{1 bar, 298 K}} + \int_{298}^T \frac{C_P}{T}~dT, \\ \int_{\textrm{1 bar}}^P V(P,T)~dP &= P V_0 \left( 1 - a + \left( a \frac{(1-b P_{th})^{1-c} - (1 + b(P-P_{th}))^{1-c}}{b (c-1) P} \right) \right)\end{split}\]

The heat capacity at one bar is given by an empirical polynomial fit to experimental data

\[C_p = a + bT + cT^{-2} + dT^{-0.5}\]

The entropy at high pressure and temperature can be calculated by differentiating the expression for \(\mathcal{G}\) with respect to temperature

\[\begin{split}\mathcal{S}(P,T) &= S_{\textrm{1 bar, T}} + \frac{\partial \int V dP }{\partial T}, \\ \frac{\partial \int V dP }{\partial T} &= V_0 \alpha_0 K_0 a \frac{C_{V0}(T)}{C_{V0}(T_\textrm{{ref}})} ((1+b(P-P_{th}))^{-c} - (1-bP_{th})^{-c} )\end{split}\]

Finally, the enthalpy at high pressure and temperature can be calculated

\[\mathcal{H}(P,T) = \mathcal{G}(P,T) + T\mathcal{S}(P,T)\]

SLB2005

The Debye model yields the Helmholtz free energy and entropy due to lattice vibrations

\[\begin{split}\mathcal{G} &= \mathcal{F} + PV, \\ \mathcal{F} &= nRT \left(3 \ln( 1 - e^{-\frac{\theta}{T}}) - \int_{0}^{\frac{\theta}{T}}\frac{\tau^3}{(e^{\tau}-1)}d\tau \right), \\ \mathcal{S} &= nR \left(4 \int_{0}^{\frac{\theta}{T}}\frac{\tau^3}{(e^{\tau}-1)}d\tau - 3 \ln(1-e^{-\frac{\theta}{T}}) \right), \\ \mathcal{H} &= \mathcal{G} + T\mathcal{S}\end{split}\]

Property modifiers

The thermodynamic models above consider the effects of strain and quasiharmonic lattice vibrations on the free energies of minerals at given temperatures and pressures. There are a number of additional processes, such as isochemical order-disorder and magnetic effects which also contribute to the total free energy of a phase. Corrections for these additional processes can be applied in a number of different ways. Burnman currently includes implementations of the following:

  • Linear excesses (useful for DQF modifications for [HollandPowell11])
  • Tricritical Landau model (two formulations)
  • Bragg-Williams model
  • Magnetic excesses

In all cases, the excess Gibbs free energy \(\mathcal{G}\) and first and second partial derivatives with respect to pressure and temperature are calculated. The thermodynamic properties of each phase are then modified in a consistent manner; specifically:

\[\begin{split}\mathcal{G} = \mathcal{G}_o + \mathcal{G}_m, \\ \mathcal{S} = \mathcal{S}_o - \frac{\partial \mathcal{G}}{\partial T}_m, \\ \mathcal{V} = \mathcal{V}_o + \frac{\partial \mathcal{G}}{\partial P}_m, \\ K_T = \mathcal{V} / \left( \frac{\mathcal{V}_o}{K_To} - \frac{\partial^2\mathcal{G}}{\partial P^2} \right)_m, \\ C_p = C_po - T \frac{\partial ^2 \mathcal{G}}{\partial T^2}_m, \\ \alpha = \left( \alpha_o \mathcal{V}_o + \frac{\partial^2 \mathcal{G}}{\partial P \partial T}_m \right) / \mathcal{V}, \\ \mathcal{H} = \mathcal{G} + T \mathcal{S}, \\ \mathcal{F} = \mathcal{G} - P \mathcal{V}, \\ C_v = C_p - \mathcal{V} T \alpha^2 K_T, \\ \gamma = \frac{\alpha K_T \mathcal{V}}{C_v}, \\ K_S = K_T \frac{C_p}{C_v}\end{split}\]

Subscripts \(_o\) and \(_m\) indicate original properties and modifiers respectively. Importantly, this allows us to stack modifications such as multiple Landau transitions in a simple and straightforward manner. In the burnman code, we add property modifiers as an attribute to each mineral as a list. For example:

from burnman.minerals import SLB_2011
stv = SLB_2011.stishovite()
stv.property_modifiers = [
['landau',
{'Tc_0': -4250.0, 'S_D': 0.012, 'V_D': 1e-09}]]
['linear',
{'delta_E': 1.e3, 'delta_S': 0., 'delta_V': 0.}]]

Each modifier is a list with two elements, first the name of the modifier type, and second a dictionary with the required parameters for that model. A list of parameters for each model is given in the following sections.

Linear excesses (linear)

A simple linear correction in pressure and temperature. Parameters are ‘delta_E’, ‘delta_S’ and ‘delta_V’.

\[\begin{split}\mathcal{G} = \Delta \mathcal{E} - T \Delta \mathcal{S} + P \Delta \mathcal{V}, \\ \frac{\partial \mathcal{G}}{\partial T} = - \Delta \mathcal{S}, \\ \frac{\partial \mathcal{G}}{\partial P} = \Delta \mathcal{V}, \\ \frac{\partial^2 \mathcal{G}}{\partial T^2} = 0, \\ \frac{\partial^2 \mathcal{G}}{\partial P^2} = 0, \\ \frac{\partial^2 \mathcal{G}}{\partial T \partial P} = 0\end{split}\]

Tricritical Landau model (landau)

Applies a tricritical Landau correction to the properties of an endmember which undergoes a displacive phase transition. These transitions are not associated with an activation energy, and therefore they occur rapidly compared with seismic wave propagation. Parameters are ‘Tc_0’, ‘S_D’ and ‘V_D’.

This correction follows [Putnis92], and is done relative to the completely ordered state (at 0 K). It therefore differs in implementation from both [SLB11] and [HollandPowell11], who compute properties relative to the completely disordered state and standard states respectively. The current implementation is preferred, as the excess entropy (and heat capacity) terms are equal to zero at 0 K.

\[Tc = Tc_0 + \frac{V_D P}{S_D}\]

If the temperature is above the critical temperature, Q (the order parameter) is equal to zero, and the Gibbs free energy is simply that of the disordered phase:

\[\begin{split}\mathcal{G}_{\textrm{dis}} = -S_D \left( \left( T - Tc \right) + \frac{Tc_0}{3} \right), \\ \frac{\partial \mathcal{G}}{\partial P}_{\textrm{dis}} = V_D, \\ \frac{\partial \mathcal{G}}{\partial T}_{\textrm{dis}} = -S_D\end{split}\]

If temperature is below the critical temperature, Q is between 0 and 1. The gibbs free energy can be described thus:

\[\begin{split}Q^2 = \sqrt{\left( 1 - \frac{T}{Tc} \right)}, \\ \mathcal{G} = S_D \left((T - Tc) Q^2 + \frac{Tc_0 Q^6}{3} \right) + \mathcal{G}_{\textrm{dis}}, \\ \frac{\partial \mathcal{G}}{\partial P} = - V_D Q^2 \left(1 + \frac{T}{2 Tc} \left(1. - \frac{Tc_0}{Tc} \right) \right) + \frac{\partial \mathcal{G}}{\partial P}_{\textrm{dis}}, \\ \frac{\partial \mathcal{G}}{\partial T} = S_D Q^2 \left(\frac{3}{2} - \frac{Tc_0}{2 Tc} \right) + \frac{\partial \mathcal{G}}{\partial T}_{\textrm{dis}}, \\ \frac{\partial^2 \mathcal{G}}{\partial P^2} = V_D^2 \frac{T}{S_D Tc^2 Q^2} \left( \frac{T}{4 Tc} \left(1. + \frac{Tc_0}{Tc} \right) + Q^4 \left(1. - \frac{Tc_0}{Tc} \right) - 1 \right), \\ \frac{\partial^2 \mathcal{G}}{\partial T^2} = - \frac{S_D}{Tc Q^2} \left(\frac{3}{4} - \frac{Tc_0}{4 Tc} \right), \\ \frac{\partial^2 \mathcal{G}}{\partial P \partial T} = \frac{V_D}{2 Tc Q^2} \left(1 + \left(\frac{T}{2 Tc} - Q^4 \right) \left(1 - \frac{Tc_0}{Tc} \right) \right)\end{split}\]

Tricritical Landau model (landau_hp)

Applies a tricritical Landau correction similar to that described above. However, this implementation follows [HollandPowell11], who compute properties relative to the standard state. Parameters are ‘P_0’, ‘T_0’, ‘Tc_0’, ‘S_D’ and ‘V_D’.

It is worth noting that the correction described by [HollandPowell11] has been incorrectly used throughout the geological literature, particularly in studies involving magnetite (which includes studies comparing oxygen fugacities to the FMQ buffer (due to an incorrect calculation of the properties of magnetite). Note that even if the implementation is correct, it still allows the order parameter Q to be greater than one, which is physically impossible.

We include this implementation in order to reproduce the dataset of [HollandPowell11]. If you are creating your own minerals, we recommend using the standard implementation.

\[Tc = Tc0 + \frac{V_D P}{S_D}\]

If the temperature is above the critical temperature, Q (the order parameter) is equal to zero. Otherwise

\[Q^2 = \sqrt{\left( \frac{Tc - T}{Tc0} \right)}\]
\[\begin{split}\mathcal{G} = Tc_0 S_D \left( Q_0^2 - \frac{Q_0 ^ 6}{3} \right) \ - S_D \left( Tc Q^2 - Tc_0 \frac{Q ^ 6}{3} \right) \ - T S_D \left( Q_0^2 - Q^2 \right) + P V_D Q_0^2, \\ \frac{\partial \mathcal{G}}{\partial P} = -V_D \left( Q^2 - Q_0^2 \right), \\ \frac{\partial \mathcal{G}}{\partial T} = S_D \left( Q^2 - Q_0^2 \right), \\\end{split}\]

The second derivatives of the Gibbs free energy are only non-zero if the order parameter exceeds zero. Then

\[\begin{split}\frac{\partial^2 \mathcal{G}}{\partial P^2} = -\frac{V_D^2}{2 S_D Tc_0 Q^2}, \\ \frac{\partial^2 \mathcal{G}}{\partial T^2} = -\frac{S_D}{2 Tc_0 Q^2}, \\ \frac{\partial^2 \mathcal{G}}{\partial P \partial T} = \frac{V_D}{2 Tc_0 Q^2}\end{split}\]

Bragg-Williams model (bragg_williams)

The Bragg-Williams model is essentially a symmetric solid solution model between endmembers, with an excess configurational entropy term predicted on the basis of the specifics of order-disorder in the mineral, multiplied by some empirical factor. Expressions for the excess Gibbs free energy can be found in [HP96]. Parameters are ‘deltaH’, ‘deltaV’, ‘Wh’, ‘Wv’, ‘n’ and ‘factor’.

Magnetic model (magnetic_chs)

This model approximates the excess energy due to magnetic ordering. It was originally described in [CHS87]. The expressions used by BurnMan can be found in [Sun91]. Parameters are ‘structural_parameter’, ‘curie_temperature’[2] (zero pressure value and pressure dependence) and ‘magnetic_moment’[2] (zero pressure value and pressure dependence).

Calculating Solid Solution Properties

Many minerals can exist over a finite region of composition space. These spaces are bounded by endmembers (which may themselves not be stable), and each individual mineral can then be described as a solid solution of those endmembers. At an atomic level, different elements substitute for one another on distinct crystallographic sites in the structure. For example, low pressure silicate garnets have two distinct sites on which mixing takes place; a dodecahedral site (of which there are three per unit cell on an eight-cation basis) and octahedral site (of which there are two per unit cell). A third tetrahedral cation site (three per unit cell) is usually assumed to be occupied solely by silicon, and therefore can be ignored in solid solution calculations. The chemical formula of many low pressure garnets exist within the solid solution:

\[\textrm{[Mg,Fe,Mn,Ca]}_3\textrm{[Al,Fe,Cr]}_2\textrm{Si}_3\textrm{O}_{12}\]

We typically calculate solid solution properties by appropriate differentiation of the Gibbs Free energy, where

\[\begin{split}\mathcal{G} = \sum_i n_i \left( \mathcal{G}_i + RT \ln \alpha_i \right)\\ \alpha_i = \gamma_i \alpha_{\textrm{ideal}, i}\end{split}\]

Implemented models

Ideal solid solutions

A solid solution is not simply a mechanical mixture of its constituent endmembers. The mixing of different elements on sites results in an excess configurational entropy

\[\mathcal{S}_{\textrm{conf}} = R \ln \prod_s (X_c^s)^{\nu}\]

where \(s\) is a site in the lattice \(M\), \(c\) are the cations mixing on site \(s\) and \(\nu\) is the number of \(s\) sites in the formula unit. Solid solutions where this configurational entropy is the only deviation from a mechanical mixture are termed ideal. From this expression, we can see that

\[\alpha_{\textrm{ideal}, i} = \prod_s (X_c^s)^{\nu}\]

Symmetric solid solutions

Many real minerals are not well approximated as ideal solid solutions. Deviations are the result of elastic and chemical interactions between ions with different physical and chemical characteristics. Regular (symmetric) solid solution models are designed to account for the simplest form of deviations from ideality, by allowing the addition of excess enthalpies, non-configurational entropies and volumes to the ideal solution model. These excess terms have the matrix form [DPWH07]

\[\mathcal{G}_{\textrm{excess}} = RT \ln \gamma = p^T W p\]

where \(p\) is a vector of molar fractions of each of the \(n\) endmembers and \(W\) is a strictly upper-triangular matrix of interaction terms between endmembers. Excesses within binary systems (\(i\)-\(j\)) have a quadratic form and a maximum of \(W_{ij}/4\) half-way between the two endmembers.

Asymmetric solid solutions

Some solid solutions exhibit asymmetric excess terms. These can be accounted for with an asymmetric solid solution [DPWH07]

\[\mathcal{G}_{\textrm{excess}} = \alpha^T p (\phi^T W \phi)\]

\(\alpha\) is a vector of “van Laar parameters” governing asymmetry in the excess properties.

\[\begin{split}\phi_i &= \frac{\alpha_i p_i}{\sum_{k=1}^{n} \alpha_k p_k}, \\ W_{ij} &= \frac{2 w_{ij}}{\alpha_i + \alpha_j} \textrm{for i<j}\end{split}\]

The \(w_{ij}\) terms are a set of interaction terms between endmembers \(i\) and \(j\). If all the \(\alpha\) terms are equal to unity, a non-zero \(w\) yields an excess with a quadratic form and a maximum of \(w/4\) half-way between the two endmembers.

Subregular solid solutions

An alternative way to create asymmetric solution models is to expand each binary term as a cubic expression [HW89]. In this case,

\[\mathcal{G}_{\textrm{excess}} = \sum_i p_i p_j^2 W_{ij} + p_j p_i^2 W_{ji}\]

Note the similarity with the symmetric solution model, the primary difference being that there are not two interaction terms for each binary.

Thermodynamic and thermoelastic properties

From the preceeding equations, we can define the thermodynamic potentials of solid solutions:

\[\begin{split}\mathcal{G}_{\textrm{SS}} &= \sum_i n_i \left( \mathcal{G}_i + RT \ln \alpha_i \right)\\ \mathcal{S}_{\textrm{SS}} &= \sum_in_i\mathcal{S}_i + \mathcal{S}_{\textrm{conf}} - \frac{\partial \mathcal{G}_{\textrm{excess}}}{\partial T} \\ \mathcal{H}_{\textrm{SS}} &= \mathcal{G}_{\textrm{SS}} + T\mathcal{S}_{\textrm{SS}}\\ V_{\textrm{SS}} &= \sum_in_iV_i + \frac{\partial \mathcal{G}_{\textrm{excess}}}{\partial P}\end{split}\]

We can also define the derivatives of volume with respect to pressure and temperature

\[\begin{split}\alpha_{P,\textrm{SS}} &= \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_P = \left( \frac{1}{V_{\textrm{SS}}}\right)\left( \sum_i\left(n_i\,\alpha_i\,V_i \right) \right) \\ K_{T,\textrm{SS}} &= V\left( \frac{\partial P}{\partial V} \right)_T = V_{\textrm{SS}} \left( \frac{1}{\sum_i\left(n_i \frac{V_{i}}{K_{Ti}} \right)} + \frac{\partial P}{\partial V_{\textrm{excess}}} \right)\end{split}\]

Making the approximation that the excess entropy has no temperature dependence

\[\begin{split}C_{P,\textrm{SS}} &= \sum_in_iC_{Pi}\\ C_{V, \textrm{SS}} &= C_{P,\textrm{SS}} - V_{\textrm{SS}}\,T\,\alpha_{\textrm{SS}}^{2}\,K_{T,\textrm{SS}} \\ K_{S,\textrm{SS}} &= K_{T,\textrm{SS}} \,\frac{C_{P,\textrm{SS}}}{C_{V,\textrm{SS}}}\\ \gamma_{\textrm{SS}} &= \frac{\alpha_{\textrm{SS}}\,K_{T,\textrm{SS}}\,V_{\textrm{SS}}}{C_{V, \textrm{SS}}}\end{split}\]

Including order-disorder

Order-disorder can be treated trivially with solid solutions. The only difference between mixing between ordered and disordered endmembers is that disordered endmembers have a non-zero configurational entropy, which must be accounted for when calculating the excess entropy within a solid solution.

Including spin transitions

The regular solid solution formalism should provide an elegant way to model spin transitions in phases such as periclase and bridgmanite. High and low spin iron can be treated as different elements, providing distinct endmembers and an excess configurational entropy. Further excess terms can be added as necessary.

Calculating Multi-phase Composite Properties

Averaging schemes

After the thermoelastic parameters (\(K_S\), \(G\), \(\rho\)) of each phase are determined at each pressure and/or temperature step, these values must be combined to determine the seismic velocity of a multiphase assemblage. We define the volume fraction of the individual minerals in an assemblage:

\[\nu_i = n_i \frac{V_i}{V},\]

where \(V_i\) and \(n_i\) are the molar volume and the molar fractions of the \(i\) th individual phase, and \(V\) is the total molar volume of the assemblage:

(11)\[V = \sum_i n_i V_i.\]

The density of the multiphase assemblage is then

(12)\[\rho = \sum_i \nu_i \rho_i = \frac{1}{V}\sum_i {n_i \mu_i},\]

where \(\rho_i\) is the density and \(\mu_i\) is the molar mass of the \(i\) th phase.

Unlike density and volume, there is no straightforward way to average the bulk and shear moduli of a multiphase rock, as it depends on the specific distribution and orientation of the constituent minerals. BurnMan allows several schemes for averaging the elastic moduli: the Voigt and Reuss bounds, the Hashin-Shtrikman bounds, the Voigt-Reuss-Hill average, and the Hashin-Shtrikman average [WDOConnell76].

The Voigt average, assuming constant strain across all phases, is defined as

(13)\[X_V = \sum_i \nu_i X_i,\]

where \(X_i\) is the bulk or shear modulus for the \(i\) th phase. The Reuss average, assuming constant stress across all phases, is defined as

(14)\[X_R = \left(\sum_i \frac{\nu_i}{X_i} \right)^{-1}.\]

The Voigt-Reuss-Hill average is the arithmetic mean of Voigt and Reuss bounds:

(15)\[X_{VRH} = \frac{1}{2} \left( X_V + X_R \right).\]

The Hashin-Shtrikman bounds make an additional assumption that the distribution of the phases is statistically isotropic and are usually much narrower than the Voigt and Reuss bounds [WDOConnell76]. This may be a poor assumption in regions of Earth with high anisotropy, such as the lowermost mantle, however these bounds are more physically motivated than the commonly-used Voigt-Reuss-Hill average. In most instances, the Voigt-Reuss-Hill average and the arithmetic mean of the Hashin-Shtrikman bounds are quite similar with the pure arithmetic mean (linear averaging) being well outside of both.

It is worth noting that each of the above bounding methods are derived from mechanical models of a linear elastic composite. It is thus only appropriate to apply them to elastic moduli, and not to other thermoelastic properties, such as wave speeds or density.

Computing seismic velocities

Once the moduli for the multiphase assemblage are computed, the compressional (\(P\)), shear (\(S\)) and bulk sound (\(\Phi\)) velocities are then result from the equations:

(16)\[V_P = \sqrt{ \frac{K_S + \frac{4}{3} G} {\rho} }, \qquad V_S = \sqrt{ \frac{G}{\rho} }, \qquad V_\Phi = \sqrt{ \frac{K_S}{\rho} }.\]

To correctly compare to observed seismic velocities one needs to correct for the frequency sensitivity of attenuation. Moduli parameters are obtained from experiments that are done at high frequencies (MHz-GHz) compared to seismic frequencies (mHz-Hz). The frequency sensitivity of attenuation causes slightly lower velocities for seismic waves than they would be for high frequency waves. In BurnMan one can correct the calculated acoustic velocity values to those for long period seismic tomography [MA81]:

\[V_{S/P}=V_{S/P}^{\mathrm{uncorr.}}\left(1-\frac{1}{2}\cot(\frac{\beta\pi}{2})\frac{1}{Q_{S/P}}(\omega)\right).\]

Similar to [MBR+07], we use a \(\beta\) value of 0.3, which falls in the range of values of \(0.2\) to \(0.4\) proposed for the lower mantle (e.g. [KS90]). The correction is implemented for \(Q\) values of PREM for the lower mantle. As \(Q_S\) is smaller than \(Q_P\), the correction is more significant for S waves. In both cases, though, the correction is minor compared to, for example, uncertainties in the temperature (corrections) and mineral physical parameters. More involved models of relaxation mechanisms can be implemented, but lead to the inclusion of more poorly constrained parameters, [MB07]. While attenuation can be ignored in many applications [TVV01], it might play a significant role in explaining strong variations in seismic velocities in the lowermost mantle [DGD+12].

User input

Mineralogical composition

A number of pre-defined minerals are included in the mineral library and users can create their own. The library includes wrapper functions to include a transition from the high-spin mineral to the low-spin mineral [LSMM13] or to combine minerals for a given iron number.

Standard minerals – The ‘standard’ mineral format includes a list of parameters given in the above table. Each mineral includes a suggested EoS with which the mineral parameters are derived. For some minerals the parameters for the thermal corrections are not yet measured or calculated, and therefore the corrections can not be applied. An occasional mineral will not have a measured or calculated shear moduli, and therefore can only be used to compute densities and bulk sound velocities. The mineral library is subdivided by citation. BurnMan includes the option to produce a LaTeX; table of the mineral parameters used. BurnMan can be easily setup to incorporate uncertainties for these parameters.

Minerals with a spin transition – A standard mineral for the high spin and low spin must be defined separately. These minerals are “wrapped,” so as to switch from the high spin to high spin mineral at a give pressure. While not realistic, for the sake of simplicity, the spin transitions are considered to be sharp at a given pressure.

Minerals depending on Fe partitioning – The wrapper function can partition iron, for example between ferropericlase, fp, and perovskite, pv. It requires the input of the iron mol fraction with regards to Mg, \(X_\mathrm{fp}\) and \(X_\mathrm{pv}\), which then defines the chemistry of an Mg-Fe solid solution according to (\(\mathrm{Mg}_{1-X_{\mathrm{Fe}}^{\mathrm{fp}}},\mathrm{Fe}_{X_{\mathrm{Fe}}^{\mathrm{fp}}})\mathrm{O}\) or \((\mathrm{Mg}_{1-X_{\mathrm{Fe}}^{\mathrm{pv}}},\mathrm{Fe}_{X_{\mathrm{Fe}}^{\mathrm{pv}}})\mathrm{SiO_3}\). The iron mol fractions can be set to be constant or varying with P and T as needed. Alternatively one can calculate the iron mol fraction from the distribution coefficient \(K_D\) defined as

(17)\[K_{D} = \frac{X_{\mathrm{Fe}}^{\mathrm{pv}}/X_{\mathrm{Mg}}^{\mathrm{pv}}}{X_{\mathrm{Fe}}^{\mathrm{fp}}/X_{\mathrm{Mg}}^{\mathrm{fp}}}.\]

We adopt the formalism of [NFR12] choosing a reference distribution coefficient \(K_{D0}\) and standard state volume change (\(\Delta \upsilon^{0}\)) for the Fe-Mg exchange between perovskite and ferropericlase

(18)\[K_{D}={K_D}_0 \:\exp\left(\frac{(P_0-P)\Delta \upsilon^{0}}{RT}\right),\]

where \(R\) is the gas constant and \(P_0\) the reference pressure. As a default, we adopt the average \(\Delta \upsilon^{0}\) of [NFR12] of \(2\cdot10^{-7}\) \(m^3 mol^{-1}\) and suggest using their \({K_D}_0\) value of \(0.5\).

The multiphase mixture of these minerals can be built by the user in three ways:

1. Molar fractions of an arbitrary number of pre-defined minerals, for example mixing standard minerals mg_perovskite (\(\mathrm{MgSiO_3}\)), fe_perovskite (\(\mathrm{FeSiO_3}\)), periclase (\(\mathrm{MgO}\)) and wüstite (\(\mathrm{FeO}\)).

2. A two-phase mixture with constant or (\(P,T\)) varying Fe partitioning using the minerals that include Fe-dependency, for example mixing \(\mathrm{(Mg,Fe)SiO_3}\) and \(\mathrm{(Mg,Fe)O}\) with a pre-defined distribution coefficient.

3. Weight percents (wt%) of (Mg, Si, Fe) and distribution coefficient (includes (P,T)-dependent Fe partitioning). This calculation assumes that each element is completely oxidized into its corresponding oxide mineral (\(\mathrm{MgO}\), \(\mathrm{FeO}\), \(\mathrm{SiO_2}\)) and then combined to form iron-bearing perovskite and ferropericlase taking into account some Fe partition coefficient.

Geotherm

Unlike the pressure, the temperature of the lower mantle is relatively unconstrained. As elsewhere, BurnMan provides a number of built-in geotherms, as well as the ability to use user-defined temperature-depth relationships. A geotherm in BurnMan is an object that returns temperature as a function of pressure. Alternatively, the user could ignore the geothermal and compute elastic velocities for a range of temperatures at any give pressure.

Currently, we include geotherms published by [BS81] and [And82]. Alternatively one can use an adiabatic gradient defined by the thermoelastic properties of a given mineralogical model. For a homogeneous material, the adiabatic temperature profile is given by integrating the ordinary differential equation (ODE)

(19)\[\left(\frac{\text{d}T}{\text{d}P}\right)_S = \frac{\gamma T}{K_S}.\]

This equation can be extended to multiphase composite using the first law of thermodynamics to arrive at

(20)\[\left(\frac{\text{d}T}{\text{d}P}\right)_S = \frac{ T \displaystyle\sum_{i} \frac{ n_i C_{Pi} \gamma_i }{K_{Si}}}{ \displaystyle\sum_{i} n_i C_{Pi} },\]

where the subscripts correspond to the \(i\) th phase, \(C_P\) is the heat capacity at constant pressure of a phase, and the other symbols are as defined above. Integrating this ODE requires a choice in anchor temperature (\(T_0\)) at the top of the lower mantle (or including this as a parameter in an inversion). As the adiabatic geotherm is dependent on the thermoelastic parameters at high pressures and temperatures, it is dependent on the equation of state used.

Seismic Models

BurnMan allows for direct visual and quantitative comparison with seismic velocity models. Various ways of plotting can be found in the examples. Quantitative misfits between two profiles include an L2-norm and a chi-squared misfit, but user defined norms can be implemented. A seismic model in BurnMan is an object that provides pressure, density, and seismic velocities (\(V_P, V_\Phi, V_S\)) as a function of depth.

To compare to seismically constrained profiles, BurnMan provides the 1D seismic velocity model PREM [DA81]. One can choose to evaluate \(V_P, V_\Phi, V_S, \rho, K_S\) and/or \(G\). The user can input their own seismic profile, an example of which is included using AK135 [KEB95].

Besides standardized 1D radial profiles, one can also compare to regionalized average profiles for the lower mantle. This option accommodates the observation that the lowermost mantle can be clustered into two regions, a ‘slow’ region, which represents the so-called Large Low Shear Velocity Provinces, and ‘fast’ region, the continuous surrounding region where slabs might subduct [LCDR12]. This clustering as well as the averaging of the 1D model occurs over five tomographic S wave velocity models (SAW24B16: [MegninR00]; HMSL-S: [HMSL08]; S362ANI: [KED08]; GyPSuM: [SFBG10]; S40RTS: [RDvHW11]). The strongest deviations from PREM occur in the lowermost 1000 km. Using the ‘fast’ and ‘slow’ S wave velocity profiles is therefore most important when interpreting the lowermost mantle. Suggestion of compositional variation between these regions comes from seismology [TRCT05][HW12] as well as geochemistry [DCT12][JCK+10]. Based on thermo-chemical convection models, [SDG11] also show that averaging profiles in thermal boundary layers may cause problems for seismic interpretation.

We additionally apply cluster analysis to and provide models for P wave velocity based on two tomographic models (MIT-P08: [LvdH08]; GyPSuM: [SMJM12]). The clustering results correlate well with the fast and slow regions for S wave velocities; this could well be due to the fact that the initial model for the P wave velocity models is scaled from S wave tomographic velocity models. Additionally, the variations in P wave velocities are a lot smaller than for S waves. For this reason using these adapted models is most important when comparing the S wave velocities.

While interpreting lateral variations of seismic velocity in terms of composition and temperature is a major goal [TDRY04][MCD+12], to determine the bulk composition the current challenge appears to be concurrently fitting absolute P and S wave velocities and incorporate the significant uncertainties in mineral physical parameters).