Skip to content

Commit

Permalink
address editorial comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@bobmyhill
bobmyhill committed Jun 22, 2023
1 parent edb4e08 commit c4aa883
Show file tree
Hide file tree
Showing 2 changed files with 123 additions and 139 deletions.
31 changes: 30 additions & 1 deletion paper/paper.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,18 @@ @article{deCapitani:2010
eprint = {https://pubs.geoscienceworld.org/msa/ammin/article-pdf/95/7/1006/3626686/14\_3354DeCapitani.pdf},
}

@article{deKoker:2013,
author = {{de Koker}, N. and {Karki}, B.~B. and {Stixrude}, L.},
title = {{Thermodynamics of the MgO-SiO$_{2}$ liquid system in Earth's lowermost mantle from first principles}},
journal = {Earth and Planetary Science Letters},
year = 2013,
month = jan,
volume = 361,
pages = {58-63},
doi = {10.1016/j.epsl.2012.11.026},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

@article{Stixrude:2022,
author = {Stixrude, Lars and Lithgow-Bertelloni, Carolina},
title = "{Thermal expansivity, heat capacity and bulk modulus of the mantle}",
Expand Down Expand Up @@ -519,4 +531,21 @@ @article{Nowak:1991
author={Nowak, U and Weimann, L},
journal={Zuse Institute Berlin, ZIB},
year={1991}
}
}

@article{Diamond:2016,
author = {{Diamond}, Steven and {Boyd}, Stephen},
title = "{CVXPY: A Python-Embedded Modeling Language for Convex Optimization}",
journal = {arXiv e-prints},
keywords = {Mathematics - Optimization and Control},
year = 2016,
month = mar,
eid = {arXiv:1603.00943},
pages = {arXiv:1603.00943},
doi = {10.48550/arXiv.1603.00943},
archivePrefix = {arXiv},
eprint = {1603.00943},
primaryClass = {math.OC},
adsurl = {https://ui.adsabs.harvard.edu/abs/2016arXiv160300943D},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
231 changes: 93 additions & 138 deletions paper/paper.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,128 +55,113 @@ bibliography: paper.bib
Many branches of physics, chemistry and Earth sciences build
complex materials from simpler constructs: rocks are composites
made up of one or more phases; phases are solutions of more than one
endmember and endmembers are usually mixtures of more than one element.
endmember; and endmembers are usually mixtures of more than one element.
The properties of the endmember building blocks at different pressures and
temperatures can be provided by a wide array of different equations of state.
The averaging of the endmember properties within solutions and composite
materials can also be achieved in several different ways. Once calculated,
temperatures rovided modelled using a wide array of different equations of state.
There are also many models for the averaging of endmember properties
within solutions and composite materials. Once calculated,
the physical properties of composite materials can be used in many
different ways.


`BurnMan` is an open source, extensible mineral physics Python module
written in Python to provide a well-tested, benchmarked toolbox that
is able to use several different methods to calculate the
physical properties of natural materials.
`BurnMan` is an open source, extensible mineral physics module
written in Python. It implements several different methods to
calculate the physical properties of natural materials.
The toolbox has a class-based, modular
design that allows users to calculate many low-level properties that
are not accessible using existing codes, and to combine various tools in
are not easily accessed using existing codes, and to combine various tools in
novel, creative ways. The module includes:

* over 10 static and thermal equations of state for endmembers
* over a dozen static and thermal equations of state for pure phases;
* commonly-used solution model formalisms (ideal, (a)symmetric, subregular)
and an option for user-defined functional forms
* popular endmember and solution datasets
(including [@Holland:2011] and [@Stixrude:2011])
and a formalism that allows users to define their own excess energy functions;
* popular endmember and solution datasets for solids and melts,
including @Holland:2011, @deKoker:2013 and @Stixrude:2022;
* an anisotropic equation of state [@Myhill:2022];
* a consistent method for combining phases into a composite assemblage,
with seismic averaging schemes including Voigt, Reuss, Voigt-Reuss-Hill
and the Hishin-Shtrikman bounds
and the Hishin-Shtrikman bounds;
* a common set of methods to output thermodynamic and thermoelastic
properties for all materials
* an equilibrium reaction solver for composites
properties for all materials;
* an solver to chemically equilibrate composite materials;
* optimal least squares fitting routines for multivariate experimental
data with (potentially correlated) errors. Examples include simultaneous
fitting of volumes, seismic velocities and enthalpies
* a "Planet" class, which self-consistently calculates gravity-pressure-density
profiles, mass, moment of inertia and temperature structure of planets given
appropriate chemical and dynamic constraints
data with (potentially correlated) errors. These allow (for example) simultaneous
fitting of pure phase and solution model parameters to experimental volumes,
seismic velocities and enthalpies of formation;
* "Planet" and "Layer" classes that self-consistently calculate
gravity, pressure, density, mass, moment of inertia and seismic velocity profiles
given chemical, thermal and dynamic constraints;
* geothermal profiles from the literature as well as the option to calculate
adiabatic profiles based on mineral assemblage
adiabatic profiles based on mineral assemblage;
* a set of high-level functions which create files readable by
seismological and geodynamic software, including: Mineos [@Masters:2011],
AxiSEM [@NissenMeyer:2014] and
ASPECT [@Kronbichler:2012;@aspect-doi-v2.4.0;@aspectmanual]
ASPECT [@Kronbichler:2012;@aspect-doi-v2.4.0;@aspectmanual]; and
* a `Composition` class, which provides a framework to convert between mass, molar,
and elemental compositions, convert to different chemical component systems,
and add or subtract components.

The project also includes a multipart tutorial
(\url{https://burnman.readthedocs.io/en/latest/tutorial.html}),
a over 40 annotated examples,
The project also includes over 40 annotated examples,
an extensive suite of unit tests and benchmarks, and
a directory containing user-contributed code from published papers.
Use of the code requires only moderate Python skills, and its modular nature
means that it can easily be customised.
The `BurnMan` project can be found at
\url{https://github.com/geodynamics/burnman}.
A multipart tutorial illustrates key functionality, including the
functions required to create the figures in this paper
(\url{https://burnman.readthedocs.io/en/latest/tutorial.html}).
Using `BurnMan` requires only moderate Python skills, and its modular nature
means that it can easily be customised.

# Statement of need
Earth Scientists are interested in a number of different material properties,
including seismic velocities, heat capacities and densities as functions of
pressure and temperature. Many of these properties are connected to each
other by physical laws. Building models of individual phases to compute
these properties can be time-consuming and prone to error,
so it is desirable to have well-tested and benchmarked software that
provides convenient functions to calculate the properties of complex
Earth, planetary and materials scientists are interested in a number
of different material properties, including seismic velocities,
densities and heat capacities as functions of pressure and temperature.
Many of these properties are connected to each
other by physical laws such as Maxwell's relations.
Building models of individual phases to compute
these properties can be time-consuming and prone to error.
It is desirable to have well-tested and benchmarked software that
makes it easy to calculate the properties of complex
composite materials from existing models, and to parameterize
new models from experimental data.
Furthermore, there are many common scientific workflows that require
physical properties as input.
thermodynamic and thermoelastic properties as input.
These are the needs satisfied by the `BurnMan` module.

# The BurnMan project
When `BurnMan` was first released [@Cottaar:2014], its focus was on
the seismic properties of the lower mantle, using a single endmember
mineral database [@Stixrude:2011] as a foundation.
This initial release had no solution model implementation,
no model for Gibbs energy or other thermodynamic potentials,
no higher level functionality for fitting of experimental data or
inversions, and no functions designed to model planetary bodies.
Since then, its scope has expanded considerably
The focus of `BurnMan` was originally on the seismic properties of the
lower mantle [@Cottaar:2014]. Its scope has now expanded
to encompass the thermodynamic and thermoelastic properties of any
geological and planetary materials
(see \url{https://github.com/geodynamics/burnman/releases}
for headline improvements). `BurnMan` now contains equations
of state for minerals and melts from several published datasets.
A common set of methods for all equations of state allows easy
access to many thermodynamic properties, including anisotropic properties
[@Myhill:2022]. A simple example is shown in \autoref{fig:qtzproperties}.
See \url{https://burnman.readthedocs.io/en/latest/tutorial_01_material_classes.html#}
for more details and the code required to make the plot.
for the history of improvements). Pure phase equations of state are designed
to be sufficiently flexible to model real-world materials from the Earth's core
to the shallowest parts of the crust (e.g. \autoref{fig:qtzproperties}).
Solution model formulations with varying levels
of non-ideality are included (e.g. \autoref{fig:garnetsolution}),
including both Gibbs and Helmholtz formulations [@Myhill:2018].
Functions are provided to convert solution models from one
endmember basis to another [@Myhill:2021]. A `Composite` class
allows calculation of the properties of assemblages
containing several phases and includes several seismic averaging schemes.

![Heat capacity and bulk sound velocities of quartz through the alpha-beta
quartz transition as found in [@Stixrude:2011]. This transition is modelled
via a Landau-type model.
\label{fig:qtzproperties}](figures/quartz_properties.png)

Several other object classes in `BurnMan` build on the endmember
equations of state. Several solution model formulations are included
(with one example model shown in \autoref{fig:garnetsolution}, see
\url{https://burnman.readthedocs.io/en/latest/tutorial_01_material_classes.html#}
for more details),
including both Gibbs and Helmholtz formulations [@Myhill:2018].
Functions are also provided to convert solution models from one
endmember basis to another [@Myhill:2021]. A composite material model
allows calculation of the properties of assemblages of several phases
and includes several seismic averaging schemes.

![Properties of pyrope-grossular garnet at 1 GPa according to a published
model [@Jennings:2015], as output by `BurnMan`.
The excess Gibbs energy is useful for calculating phase equilibria by
Gibbs minimization, while the endmember activities can be used to
determine equilibrium via the equilibrium relations [@Holland:1998].
\label{fig:garnetsolution}](figures/mg_ca_gt_properties.png)

Another class implemented in `BurnMan` is the `Composition` class
which provides a framework to convert between mass, molar,
and elemental compositions, convert to different component systems,
and add or subtract chemical components (see
\url{https://burnman.readthedocs.io/en/latest/tutorial_02_composition_class.html#}).

`BurnMan` also includes planetary `Layer` and `Planet` classes,
which can be used to construct planetary models with self-consistent
`BurnMan` also includes planetary `Layer` and `Planet` classes
that can be used to construct planetary models with self-consistent
pressure, gravity and density profiles
and calculate seismic properties through those bodies
\autoref{fig:zog} shows a 1D profile through Planet Zog, a planet much like Earth
(see \url{https://burnman.readthedocs.io/en/latest/tutorial_03_layers_and_planets.html}
for more details).
Tools are provided to compare those seismic properties with published seismic models of
the Earth, and to produce input files to compute synthetic
and calculate seismic properties through those bodies.
\autoref{fig:zog} shows the output from a model that resembles Earth.
Tools are provided to compare predicted seismic properties
with published seismic models of the Earth,
and to produce input files to compute synthetic
seismic data using other codes, including
AxiSEM [@NissenMeyer:2014] and Mineos [@Woodhouse:1988;@Masters:2011].

Expand All @@ -199,25 +184,15 @@ The computed geotherm is compared to several from the literature
\label{fig:zog}](figures/zog.png)


Separate from the core classes within `BurnMan`, the module also
includes many utility functions. These include functions to
fit thermodynamic model parameters for endmembers
and solutions including full error propagation
(\autoref{fig:fit}, \url{https://burnman.readthedocs.io/en/latest/tutorial_04_fitting.html#}).


![Optimized fit of a PV equation of state [@Holland:2011] to
stishovite data [@Andrault:2003], including 95% confidence intervals
on both the volume and the bulk modulus.
\label{fig:fit}](figures/stishovite_fitting.png)


The utility functions `fit_composition_to_solution()` and
`fit_phase_proportions_to_bulk_composition()`
use weighted constrained least squares to estimate endmember or
phase proportions and their corresponding covariances, given a bulk composition.
The fitting functions apply appropriate
non-negativity constraints
`BurnMan` also includes many utility functions. These include functions
that fit parameters for pure phase
and solution models to experimental data
including full error propagation (\autoref{fig:fit}).
Other fitting functions include `fit_composition_to_solution()` and
`fit_phase_proportions_to_bulk_composition()` that
use weighted constrained least squares using cvxpy [@Diamond:2016] to
estimate endmember or phase proportions given a bulk composition.
These fitting functions apply appropriate non-negativity constraints
(i.e. that no species can have negative proportions on a site,
and that no phase can have a negative abundance
in the bulk). An example of `fit_phase_proportions_to_bulk_composition()`
Expand All @@ -226,15 +201,11 @@ Loss of an independent component from the bulk composition
can be tested by adding another phase with the composition of
that component (e.g. Fe) and checking that the amount of that phase
is zero within uncertainties.
Phase proportion covariances $C$ are given as
a function of the bulk composition covariances $K$,
and the phase compositions $A$
$$C_{im} = (A_{ji} K^{-1}_{jk} A_{km})^{-1}.$$
The weighted residuals are calculated
using the phase proportions $x$ and the bulk composition $b$:
$$\chi^2 = (K^{-0.5}_{ij} A_{jk} x_k - K^{-0.5}_{ik} b_k)(K^{-0.5}_{ip} A_{pq} x_q - K^{-0.5}_{iq} b_q).$$
See \url{https://burnman.readthedocs.io/en/latest/tutorial_04_fitting.html#}
for more details and the code required to make \autoref{fig:mars_fit}.

![Optimized fit of a PV equation of state [@Holland:2011] to
stishovite data [@Andrault:2003], including 95% confidence intervals
on both the volume and the bulk modulus.
\label{fig:fit}](figures/stishovite_fitting.png)

![Mineral phase proportions in the mantle of Mars, estimated by using
the method of constrained least squares on high pressure experimental
Expand All @@ -252,46 +223,30 @@ HeFESTo [@Stixrude:2022] and
FactSAGE [@Bale:2002].
Instead, it provides two methods to deal with the problem of
thermodynamic equilibrium: (1) reading in a pressure-temperature table
of precalculated properties into a Material class
(from which derivative properties can be calculated),
and (2) an function called `equilibrate()`
that chemically equilibrates a known assemblage under equality constraints.
The equilibrate function is similar to that implemented in THERMOCALC
[@Holland:1998]. In this function, the user chooses an assemblage
of precalculated properties into a Material class that allows derivative
properties to be calculated, and (2) a function called `equilibrate()`
that equilibrates a known assemblage under user-defined constraints.
This function requires an assemblage
(e.g. olivine, garnet and orthopyroxene), a starting bulk composition,
desired equality constraints, and optionally one or more
compositional degrees of freedom. The equilibrate function attempts
to find the remaining unknowns that satisfy those constraints using a
damped Newton root finder [@Nowak:1991], modified to allow constraints
using the method of Lagrange multipliers.

There are a number of equality constraints implemented in `BurnMan`.
These include fixed pressure, temperature, entropy or volume,
"PT_ellipse" (which finds an equilibrium point on an ellipse with center given
by a fixed pressure and temperature), a fixed molar fraction
of a particular phase (which could be zero), or a compositional constraint on a
solution phase (such as a certain ratio of Mg and Fe on a particular site).
compositional degrees of freedom. The equilibrate function solves the
equilibrium relations [@Holland:1998] using a
damped Newton root finder [@Nowak:1991].

The `equilibrate()` function allows the user to select from a number of equality
constraints, including fixed pressure, temperature, entropy or volume, or
compositional equalities such as a fixed molar fraction of a phase,
or a certain ratio of Mg and Fe on a particular site.
The number of constraints required is two at fixed bulk composition,
and one more for each degree of compositional freedom.

An example of use of the equilibrate function is shown in \autoref{fig:eqm}.
See \url{https://burnman.readthedocs.io/en/latest/tutorial_05_equilibrium.html#}
for more details and the code required to make \autoref{fig:eqm}. The manual
also contains significantly more detail on this function.
and one more for each degree of compositional freedom. An example of the use
of the equilibrate function is shown in \autoref{fig:eqm}.
Full details may be found in the manual and tutorial.

![The olivine phase diagram at three different temperatures as computed
using the equilibrate routines in `BurnMan`. The solution model properties
are taken from the literature [@Stixrude:2011].
\label{fig:eqm}](figures/olivine_equilibrium.png)

The features used in `BurnMan` are documented in the codebase
(\url{https://github.com/geodynamics/burnman}), and in the
manual (\url{https://burnman.readthedocs.io}). A Jupyter notebook tutorial
(available at \url{https://geodynamics.github.io/burnman/#try})
introduces some of the core functionality, including demonstrations that
recreate the figures presented here. More features are investigated
in the extensive set of examples bundled with the source code.

# Past and ongoing research projects
In addition to mantle studies
[@Cottaar:2014b;@Ballmer:2017;@Jenkins:2017;@Thomson:2019;@Houser:2020],
Expand Down

0 comments on commit c4aa883

Please sign in to comment.