A Computational Framework for Mineralogical Thermodynamics of the Earth's Mantle

Thomas C. Chust
<thomas.chust@uni-bayreuth.de>

Material Properties

Temperature field of a convection simulation from Nakagawa et al. [2010]

$$\begin{align*} \nabla\cdot \left(\color{red}{\rho} u\right) &= 0, \\ -\nabla p + \color{red}{\eta} \Delta u + \color{red}{\rho} g &= 0, \\ \color{red}{\rho} T \, \left(\frac{\partial}{\partial {t}} + u\cdot\nabla\right) \color{red}{S} &= \color{red}{k} \, \Delta T + \color{red}{\bar\tau} \, \nabla u + \color{red}{\rho H^*}. \end{align*}$$

Global seismic wave propagation simulation from Schuberth et al. [2012]

$$ \partial_t^2 x = \color{red}{v}^2 \Delta x, \\ \color{red}{v_p} = \sqrt{\frac{\color{red}{\kappa_S} + \frac43\color{red}{\mu}}{\color{red}{\rho}}}, \color{red}{v_s} = \sqrt{\frac{\color{red}{\mu}}{\color{red}{\rho}}}. $$

Both simulations and data interpretation in the geosciences need information about material properties of mantle minerals.

For example, you can see the temperature field from a mantle convection simulation on the left POINT.

The governing equations for this problem describe POINT mass, POINT momentum and POINT energy conservation and include the density (POINT $\rho$) in several places, notably as the key ingredient in the POINT buoyancy term.

Entropy (POINT $S$) is similarly important for the energy equation to describe the heat exchange at phase transitions. Thermal conductivity (POINT $k$) describes diffusion of temperature while viscosity (POINT $\eta$) describes diffusion of mechanical tension.

On the right you can see a temporal snapshot of the wave field from a simulated earthquake POINT.

The wave equation describing this problem includes a velocity (POINT $v$) that in turn depends on elastic moduli (POINT $\kappa_s, \mu$).

Neither these properties nor the composition of the mantle can easily be measured in-situ, thus it is desirable to have a theoretical framework that allows the extrapolation of experimental knowledge to a wide range of pressure and temperature conditions and can predict both phase stability and thermoelastic properties.

Thermodynamic State

Category	Properties	Dependencies
Thermal	$S, C_p$	$\leftarrow \partial_T^{1,2} G$
Elastic	$V, \rho, \kappa_S, \color{#88f}{\mu}$	$\leftarrow \partial_{p,T}^{1,2} G$
Transport	$k, \eta$	Not linked to $G$
Potentials	$E, A, H, G$	$\mathrm{d}G = V\,\mathrm{d}p - S\,\mathrm{d}T$

The material properties we saw on the last slide can be broadly classified into three categories: POINT thermal, elastic and transport properties.

In addition, thermodynamics knows the very useful concept of potential functions. Once we know the functional form of the internal energy (POINT $E$), Helmholtz energy (POINT $A$), enthalpy (POINT $H$) or Gibbs energy (POINT $G$), we can extract the thermal and elastic properties using partial derivatives.

Note that transport properties cannot be derived from equilibrium thermodynamics. The shear modulus can be computed in a way consistent with the other thermoelastic properties, but requires additional model parameters.

The potential functions differ by their natural free variables, Gibbs energy is tied to pressure and temperature and is also particularly useful due to its connection with chemical equilibrium. To obtain a potential function, its total differential may be integrated.

Birch-Murnaghan–Mie-Grüneisen-Debye Model

Here you can see an examplary flow diagram for the computation of Gibbs energy (POINT $G(p, T)$) in a physical model popular for its good $p, T$ extrapolation behaviour.

The database parameters of the model are shown in green ellipses.

The model is formulated in terms of Helmholtz energy, with free variables temperature and volume, or rather isotropic finite strain (POINT $V(f)$), which is directly tied to volume.

To obtain Helmholtz energy, integration follows a path of

isothermal compression (POINT $\Delta f A(f, T)$), followed by
isochoric heating (POINT $\Delta T A(f, T)$).

The Helmholtz energy of the reference state (POINT $A_0$) is used as an integration constant and a Legendre transform adds $pV$ yielding Gibbs energy.

To obtain the isotropic finite strain parameter of the target state given $p, T$, the combined elastic and thermal pressure model is inverted numerically (POINT $f(p, T)$).

The elastic part of the model uses a Birch-Murnaghan equation of state to describe the dependency of volume (POINT $V(f)$) and elastic pressure (POINT $p_{\mathrm{el}}(f)$) on isotropic finite strain with the model parameters volume (POINT $V_0$), isothermal bulk modulus (POINT $K_0$) and pressure derivative of the isothermal bulk modulus (POINT $K_{0,p}$) in the reference state.

The thermal part of the model uses the Mie-Grüneisen-Debye statistical description of lattice vibrations to link thermal pressure (POINT $p_{\mathrm{th}}(f, T)$) to the database parameters Debye temperature (POINT $\theta_0$), Grüneisen parameter (POINT $\gamma_0$) and logarithmic volume derivative of the Grüneisen parameter (POINT $q_0$) in the reference state.

Programming Interface

Programming Interface

The structure of thermodynamic theory has a useful trait: We can reason about the thermodynamic potential, the thermoelastic properties of a system and their relations without a concrete physical model in mind.

Our simumation software preserves this trait in its object oriented design. POINT Two alternative implementations of concrete equations of state are provided. They remain interchangeable and are accessed through the same abstract interfaces.

Details of those other models are omitted for brevity here. Computational results presented throughout this talk use the Birch-Murnaghan–Mie-Grüneisen-Debye model together with the dataset published by Stixrude and Lithgow-Bertelloni [2011].

Programming Interface

Transitions of Perovskite Endmembers

Energetic equilibrium between the $\ce{Mg}$- (blue), $\ce{Fe}$- (red) and $\ce{Al}$-endmembers (green) of bridgmanite and post-perovskite in $p, T$ space.

Let us first look at pairs of fixed-composition phases. While other representations are possible, the database by Stixrude and Lithgow-Bertelloni [2011] contains $\ce{MgSiO3}$, $\ce{FeSiO3}$ and $\ce{Al2O3}$ endmembers for the bridgmanite and post-peroskite solutions. This plot shows lines of equal Gibbs energy, that is the phase boundaries, between matching pairs of these endmember phases. (POINT blue, red and green respectively).

The shape of the phase boundaries is radically different ranging from a steep positive clapeyron slope for the $\ce{Mg}$-endmembers to a steep negative slope for the $\ce{Al}$-endmembers. The equilibrium between the $\ce{Fe}$-endmembers is far less sensitive to temperature, but is also located at much higher pressure than for the $\ce{Mg}$-endmembers.

Significant variations in the bridgmanite–post-perovskite phase boundary depending on the $\ce{Fe}$ and $\ce{Al}$ content can be expected based on these results.

Transitions of Perovskite Endmembers

Energetic equilibrium between the $\ce{Mg}$- (blue), $\ce{Fe}$- (red) and $\ce{Al}$-endmembers (green) of bridgmanite and post-perovskite in $p, T$ space.

Chemical Equilibrium

In order to compute chemical equilibrium between multiple phases with varying composition, we need three more things:

The Gibbs energy of phases modelled as solid solutions has to be computed.
The Gibbs energy of aggregates of phases has to be computed.
An energetically optimal collection of phases has to be chosen at given $p, T$ and bulk composition.

READ SLIDE; AFTER (1): Since thermodynamic potentials are additive, a linear combination of constituent energies can be used for a solution. However, the sum is augmented by the effects of statistical mixing entropy following the Bragg-Williams approximation as well as by non-ideal binary interaction energy terms following the formulation by Hillert and Staffansson [1970].

AFTER (2): For an aggregate, a linear combination of constituent energies is sufficient.

AFTER (3): Fortunately, finding the minimum of a linear objective function like the aggregate Gibbs energy is a well understood problem that can be solved efficiently using certain linear algebra libraries.

FMS System

Phase diagram computed for $\ce{FeO-MgO-SiO2}$ composition with $\frac{\ce{Fe}}{\ce{Fe + Mg}} \approx 0.11$ as a function of $p$ and $T$.

So let's look at some phase diagrams for the Earth's mantle. The following slides use bulk compositions based on pyrolite, a uniform chemical model for the Earth's mantle that approximates its mineralogical state, rheology and dynamics. The phase assemblages are governed by the coexistence of $\ce{Mg2SiO4}$- and $\ce{MgSiO3}$-based minerals at upper mantle conditions, replaced by perovskite- and oxide-dominated assemblages in the lower mantle ($x(\mathrm{SiO_2}) \sim 0.4$).

We explored the sensitivity of phase relations to chemical composition by sequentially increasing the number of components from the binary MS ($\ce{MgO-SiO2}$) to the six-component NCFMAS ($\ce{Na2O-CaO-FeO-MgO-Al2O3-SiO2}$) system.

An isentrope with a potential temperature of $1600\,\mathrm{K}$ is displayed as a proxy for the geotherm POINT.

FMS System

Phase diagram computed for $\ce{FeO-MgO-SiO2}$ composition with $\frac{\ce{Fe}}{\ce{Fe + Mg}} \approx 0.11$ as a function of $p$ and $T$.

FMS System

Phase diagram computed for $\ce{FeO-MgO-SiO2}$ composition with $\frac{\ce{Fe}}{\ce{Fe + Mg}} \approx 0.11$ as a function of $p$ and $T$.

For the pyroxene-based assemblages, the phase diagram implies that, along the isentrope, orthopyroxene transforms to high-pressure clinopyroxene POINT. After intersecting the majorite stability field POINT, high-pressure clinopyroxene dissociates to ringwoodite and stishovite POINT before recombining into bridgmanite POINT.

Akimotoite occurs only at temperatures below the isentrope POINT. The small akimotoite stability field and large stishovite stability field in this model are not supported by experiments, e.g. Ishii et al. [2011], Ohtani et al. [1991], Ito and Yamada [1982].

At higher $T$, majorite garnet transforms to bridgmanite directly with a positive Clapeyron slope POINT.

Additional Components

With the addition of more chemical components,

new phases and their stability fields occur,
new endmembers are introduced and form mineral solutions whose stability fields generally expand, and/or
the variance of mineral assemblages increases, resulting in binary or higher-order phase loops which allow the coexistence of phases previously separated by sharp boundaries.

New Phases

Relative atomic amounts of stable phases for $\ce{FeO-MgO-SiO2}$ composition along the isentrope with a potential temperature of $1600\,\mathrm{K}$.

New Phases

Relative atomic amounts of stable phases for $\ce{CaO-FeO-MgO-SiO2}$ composition along the isentrope with a potential temperature of $1600\,\mathrm{K}$.

Stability Fields

Relative atomic amounts of stable phases for $\ce{FeO-MgO-SiO2}$ composition along the isentrope with a potential temperature of $1600\,\mathrm{K}$.

Stability Fields

Relative atomic amounts of stable phases for $\ce{FeO-MgO-Al2O3-SiO2}$ composition along the isentrope with a potential temperature of $1600\,\mathrm{K}$.

Coexistence Regions

Phase diagram for $\ce{Na2O-CaO-FeO-MgO-Al2O3-SiO2}$ composition as a function of $p$ and $T$. Lines indicate phase boundaries for $\ce{MgSiO3}$ (blue), $\ce{FeSiO3}$ (red) and $\ce{Al2O3}$ endmembers (green) of bridgmanite and post-perovskite.

As an example of emerging phase coexistence regions, consider the phase boundary between bridgmanite and post-perovskite in the lower mantle:

In the MS system, the phase boundary would match the POINT blue line indicating equilibrium between the $\ce{MgSiO3}$ endmembers of the two solid solutions.

But in the full NCFMAS system, the equilibrium between the solutions looks different:

The higher transition pressure in the $\ce{FeSiO3}$ endmembers (POINT red line) leads to a POINT coexistence region at low temperatures.

The inverse Clapeyron slope for the transition between the $\ce{Al2O3}$ endmembers (POINT green line) leads to a POINT coexistence region at high temperatures, extending in the other direction along the $p$ axis.

Mantle Adiabats

$$ \partial_pT\mid_S = \frac{\alpha V}{C_p} T $$

Self-consistently computed isentropes for two different pyrolite-based compositions in comparison with the adiabat from Katsura et al. [2010].

The phase diagram cuts on the last slides already used adiabatic temperature profiles implicitly. This plot compares adiabats with potential temperatures of $1600\,\mathrm{K}$ in the pressure range of $0-30\,\mathrm{GPa}$ explicitly.

Within a specific phase stability field, temperature along an isentrope increases smoothly due to self-compression, but at phase transitions, temperature can change discontinuously through volume collapse and latent heat release or consumption.

Our simulation software computes adiabats by actually tracing a curve of constant entropy starting from some specified point in $p, T$ space and examining the stable phase assemblages along the path.

Katsura et al. [2010] have published an adiabat based on thermoelastic properties of major mantle minerals in the FMS system and their observed phase transformations. This adiabat resembles our result for the FMS system, but our simulation encounters additional phase transitions.

However, the magnitude of temperature changes becomes smaller in the multiphase assemblages of the NCFMAS system due to competing energetic effects of transitions happening in parallel.

Conclusions

High accuracy of thermo-elastic parameters or good $p, T$ extrapolation behaviour is possible.
Calibration of phase equilibria may be specific to certain composition ranges.
Complex phase transitions buffer latent heat effects and result in smoother adiabats.
Elastic property profiles computed for homogeneous chemical composition or mechanical mixtures are very similar.

MMA-EoS

https://bitbucket.org/chust/eos

Versatile simulation software for minerals and mantle aggregates.
Computes thermoelastic properties and chemical equilibrium states.
Supports multiple widely used equations of state.
Features an easy-to-use solution model.
Offers a free, open and extensible implementation.