BurnMan – a Python toolkit for planetary geophysics, geochemistry and thermodynamics

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. The properties of the endmember building blocks at different pressures and temperatures can be 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.

• over a dozen static and thermal equations of state for pure phases; • commonly-used solution model formalisms (ideal, (a)symmetric, subregular) and a formalism that allows users to define their own excess energy functions; • popular endmember and solution datasets for solids and melts, including Holland & Powell (2011), de Koker et al. (2013 and Stixrude & Lithgow-Bertelloni (2021); • 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 Hashin-Shtrikman bounds; • a common set of methods to output thermodynamic and thermoelastic properties for all materials; • a solver to chemically equilibrate composite materials; • optimal least squares fitting routines for multivariate experimental 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; • a set of high-level functions which create files readable by seismological and geodynamic software, including: Mineos (Masters et al., 2011), AxiSEM (Nissen-Meyer et al., 2014 and ASPECT (Bangerth et al., 2022a(Bangerth et al., , 2022bKronbichler et al., 2012); 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 includes over 40 annotated examples, an extensive suite of unit tests and benchmarks, and a directory containing user-contributed code from published papers. A multipart tutorial illustrates key functionality, including the functions required to create the figures in this paper (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, 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 thermodynamic and thermoelastic properties as input. These are the needs satisfied by the BurnMan module.

The BurnMan project
The focus of BurnMan was originally on the seismic properties of the lower mantle (Cottaar, Heister, et al., 2014). Its scope has now expanded to encompass the thermodynamic and thermoelastic properties of any geological and planetary materials (see https://github.com/ geodynamics/burnman/releases 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. Figure 1). Solution model formulations with varying levels of non-ideality are included (e.g. Figure 2), including both Gibbs and Helmholtz formulations (Myhill, 2018). Functions are provided to convert solution models from one endmember basis to another (Myhill & Connolly, 2021). A Composite class allows calculation of the properties of assemblages containing several phases and includes several seismic averaging schemes.

Figure 1:
Heat capacity and bulk sound velocities of quartz through the alpha-beta quartz transition as found in (Stixrude & Lithgow-Bertelloni, 2011). This transition is modelled via a Landau-type model.

Figure 2:
Properties of pyrope-grossular garnet at 1 GPa according to a published model (Jennings & Holland, 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 & Powell, 1998).
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. Figure 3 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 (Nissen-Meyer et al., 2014) and Mineos (Masters et al., 2011;Woodhouse, 1988).

Figure 3:
A 1D profile through Planet Zog, a planet much like Earth, with an inner and outer core (blue and orange layers), isentropic convecting lower and upper mantle (green and red), and a depleted lithosphere (lilac) split into mantle and crust. The mineralogy/composition of each layer is chosen by the user. Zog has the same mass (5.972e+24 kg) and moment of inertia factor (0.3307) as Earth. BurnMan ensures that the gravity and pressure profiles satisfy hydrostatic equilibrium, and allows different layers to have different thermal profiles, including an isentropic profile with thermal boundary layers (shown here for the upper mantle, lower mantle and for the core). Depth dependent changes to density, gravity, pressure (solid blue lines) are compared with the Preliminary Reference Earth Model (PREM; dotted orange line, (Dziewonski & Anderson, 1981)). The computed geotherm is compared to several from the literature (Alfe et al., 2007;Anderson, 1982;Anzellini et al., 2013;Brown & Shankland, 1981;Stacey, 1977).
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 (Figure 4). Other fitting functions include fit_composition_to_solution() and fit_phase_proportions_to_bulk_composition() that use weighted constrained least squares using cvxpy (Diamond & Boyd, 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() that uses real experimental data (Bertka & Fei, 1997) is shown in Figure 5. 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.  (Holland & Powell, 2011) to stishovite data (Andrault et al., 2003), including 95% confidence intervals on both the volume and the bulk modulus.  (Bertka & Fei, 1997). Weighted residuals (misfits) are also shown, indicating that the majority of experimental run products are consistent with the reported starting composition.
BurnMan does not attempt to replicate Gibbs minimization codes, of which there are many, such as PerpleX (Connolly, 2009), MELTS (Ghiorso & Sack, 1995), MageMin (Riel et al., 2022), TheriakDomino (Capitani & Petrakakis, 2010), HeFESTo (Stixrude & Lithgow-Bertelloni, 2021) and FactSAGE (Bale et al., 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 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 solves the equilibrium relations (Holland & Powell, 1998) using a damped Newton root finder (Nowak & Weimann, 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 the use of the equilibrate function is shown in Figure 6. Full details may be found in the manual and tutorial.