pyMelt : An extensible Python engine for mantle melting calculations

Modelling the melting of Earth’s mantle is crucial for understanding the distribution of volcanic activity on Earth and for testing models of mantle convection and mantle lithological heterogeneity. pyMelt is a new open-source Python library for calculating the melting behaviour of multi-lithology mantle and can be used to predict a number of geophysical and petrological observations, including melt productivity, spreading centre crustal thickness, lava trace element concentrations, and olivine crystallisation temperatures. The library is designed to be easily extensible, allowing melting models to be added, different methods for calculating lava chemistry to be applied, and new melting dynamics and properties to be incorporated.


INTRODUCTION
Models for melting of Earth's convecting mantle can provide quantitative constraints on its thermal variability [e.g.Ball et al. 2021], compositional variability [e.g. Brown and Lesher 2014;Gleeson et al. 2021], and primary melt compositions [e.g.Jennings et al. 2016].These melting models calculate mantle melting behaviour either by minimising thermodynamic potentials at each calculation step [e.g.Smith and Asimow 2005], or by using expressions parameterised directly from melting experiments [e.g.Lambart et al. 2016;Krein et al. 2020].The parameterised approach is particularly useful for calculations requiring many runs of a melting model, for example when inverting for mantle properties from geochemical or geophysical observations [e.g.McKenzie and O'Nions 1991;Matthews et al. 2021].
Here, we present and describe pyMelt, an open-source extensible Python library that employs the parameterised approach, providing a powerful and flexible tool for calculating the melting behaviour of lithologically heterogeneous mantle.This library allows users to calculate melt productivity during isobaric melting or adiabatic decompression melting at an assigned mantle potential temperature   .Lithologies and melting parameterisations employed to calculate melt productivity can be user-defined or chosen from a list of published models (Section 3).In addition, calculated melt pathways can be used to estimate igneous crustal thickness, lava trace element concentrations, and olivine crystallisation temperatures at spreading centres or intra-plate settings (Section 6 and 7).More complex melting scenarios and other geochemical and geophysical predictions can be built from the existing library.
A number of software tools are available that can perform similar, though in some respects more limited, calculations, including INVMEL [McKenzie and O'Nions 1991], REEBOX-PRO [Brown and Lesher 2016], MELT-PX [Lambart et al. 2016], Petrogen [Krein et al. 2020], and BDD21 [Ball et al. 2022].While these packages have been used extensively in studies of mantle melting, pyMelt offers a number of advantages.pyMelt is the only package that is simultaneously open-source, incorporates mantle lithological heterogeneity, and can be easily integrated with other Python libraries (e.g.Monte Carlo inversion tools).
Open-source software is an essential part of open, transparent, and reproducible science, and it provides the basis for the development of more advanced codes and integration with other libraries.The importance of modelling the effects of lithological heterogeneity on mantle melting behaviour for accurately predicting magmatic productivity, compositions, and temperatures is becoming increasingly clear [Phipps Morgan 2001;Pertermann and Hirschmann 2003;Ito and Mahoney 2005;Shorttle et al. 2014;Matthews et al. 2021].pyMelt therefore occupies an important niche in mantle melting calculations and will provide the basis for solving new problems related to melt generation and magmatism for many years to come.
In this manuscript we review the main features of pyMelt, its library structure, its underlying mathematical formulation, and the computational approaches it takes.Users are directed to the pyMelt documentation † for a comprehensive guide to using the library, and to the interactive tutorials available online via myBinder ‡ .The pyMelt repository is hosted on GitHub § , where the code can be obtained, bugs reported, new features requested, and new contributions made.pyMelt can also be installed directly using the pip package manager.

Terminology
Throughout this manuscript we use Python terminology in our description of the library's implementation.Here we give a short description of these terms, but it is not necessary to understand this terminology to use pyMelt: • Class: A collection of methods and properties which serves as a template for creating an object.For example, the properties of the mantle are contained in a class.
• Method: An algorithm that belongs to a class that may use user-input, properties of the class, or other methods of the class to produce a result.For example, the solidus temperature of the mantle is calculated using a method of the mantle class.
• Object: An instance of a class which is used by the library to perform calculations.For example, pyMelt calculations are performed with a mantle object, created from the mantle class.
More extensive definitions can be found in the Python documentation.

pyMelt STRUCTURE
Figure 1 summarises the modular structure of the pyMelt library and its workflow.The melting behaviours of individual lithologies (e.g.lherzolite or pyroxenite) are contained within lithology objects (Section 3), which can be combined in specified mass fractions φ to form a Mantle object (Section 5).Any lithology object in pyMelt can be converted to a hydrousLithology object where the effect of water on its solidus and melt productivity is estimated (Section 4).An adiabatic decompression calculation can then be performed on a Mantle object at a specified mantle potential temperature   (the temperature a parcel of mantle would have following decompression to 0 GPa while undergoing no chemical changes) using its adiabaticMelt method (Section 5) which returns a meltingColumn object.Alternatively, an isobaric melting calculation can be performed at a specified temperature using the isobaricMelt method.If desired, the trace element concentrations in these melts can be calculated by calling the calculateChemistry method of the meltingColumn object (Section 6).To calculate geological setting specific properties, such as crustal thickness at a mid-ocean ridge, a geoSetting object can be created using the meltingColumn object (Section 7).The results can then be extracted from the geoSetting object and plotted, used in further calculations, or saved (see the tutorial notebooks).

MANTLE LITHOLOGIES
The experiments that are used to parameterise melting models are performed on particular bulk compositions, or lithologies.This means that each melting parameterisation, unless parameterised also for bulk composition [e.g.Lambart et al. 2016], represents a particular lithology, with its own melting behaviour.At a minimum, a lithology object has methods defined for the solidus temperature  solidus () (TSolidus), liquidus temperature  liquidus () (TLiquidus), and for the melt fraction  (, ,  solidus ,  liquidus ) (F), where  and  are pressure and temperature.The models contained within pyMelt already have methods for (/)  (dTdF) and (/)  (dTdP), but pyMelt can also calculate their values numerically from the TSolidus, TLiquidus, and F methods.Additionally, each lithology object has properties associated with it for the density of the solid lithology ρ  and its melt ρ  , their thermal expansivities α  and α  respectively, the heat capacity   , and the entropy change on melting Δ.
Table 1 lists the melting parameterisations for which lithology classes are defined in pyMelt.The pyroxenitic lithologies are categorised as having a relative excess or

Lithology classes
T Solidus (P), T Liquidus (P), ∂T/∂P(P,T), ∂T/∂F(P,T), ∆Smelt, cp, ρ s , ρ l , α s , α l  deficit in silica, with the silica-excess parameterisations in pyMelt representing a mid-ocean ridge basalt like composition, and the silica-deficient pyroxenites representing a mixture of basalt and lherzolite.The lithologies included already in pyMelt (with the exception of McKenzie and Bickle [1988]) have simple analytical expressions for  (, ); however, more complex models could be included, providing a method for  (, ) can be created.For example, the widely used parameterisations in the Petrogen software [Krein et al. 2020] could be added to pyMelt in the future, providing a method is created for numerically solving the more complex algorithms they use; however, calculations using such a model would likely be much slower to compute.

HYDROUS MELTING
Any lithology in pyMelt can be turned into a hydrous lithology using the hydrousLithology class.To approximate the effect of hydrous melting a similar formulation to that developed by Katz et al. [2003] is used.In this formulation the solidus temperature is depressed according to their Equation 16:

Lherzolite
where  and γ are constants and  H 2 O is the water concentration in the melt in wt.%.This differs slightly from the equations developed by Katz et al. [2003] as they did not apply the −  γ H 2 O term to every instance of  solidus in their expressions (e.g. the denominator of their Equation 19).The amount of water that can be dissolved in magmas before a H 2 O vapour phase is exsolved is limited, but increases with pressure.Katz et al. [2003] model this effect using: where χ 1 , χ 2 , and λ are constants.This equation is also implemented in pyMelt.
The concentration of H 2 O present in the melt decreases as melting proceeds, owing to H 2 O partitioning favourably into the melt and being continually diluted by new additions of magma.Katz et al. [2003] modelled this change by using the batch melting equation: where 1(which is then used to calculate the values of the other melting functions).
We extend this formulation to consider also the removal of H 2 O from the system by near-fractional melting, as modelled previously by Asimow et al. [2003]; however, since the melting models in pyMelt are implicitly expressions of batch (or equilibrium) melting, modelling the effect of H 2 O extraction by fractional melting cannot be done entirely self-consistently.This effect is approximated in pyMelt by replacing Equation 3with an expression for near-fractional melting: where Φ is the porosity during melting.When a hydrousLithology object is created, the supporting methods from the original lithology object are copied, along with the TLiquidus method (which is not changed by the hydrous melting extension).A new method for TSolidus is defined, which applies Equation 1 to the original TSolidus method.Since the value of  depends on  H 2 O , which itself depends on , a new F method is created, which solves the equation: where  calc () is the original F method (which will provide the hydrous melt fraction as it calls the modified TSolidus method) and  guess is the value changed by the root finding method.pyMelt uses the brentq algorithm implemented in the SciPy.optimize.root_scalarmethod [Virtanen et al. 2020].
Hydrous melting also changes the other melting functions, so new methods for and    are created, which calculate the values numerically using SciPy.misc.differentiate,alongside SciPy.optimize.root_scalarto find the - curve at constant .
The default values for , γ, χ 1 , χ 2 , λ, and  H 2 O are taken from Katz et al. [2003] who calibrated them for hydrouslherzolite melting.While pyMelt provides the opportunity to model hydrous-pyroxenite melting in the same way, the user must choose appropriate constant values for pyroxenite.

MANTLE MELTING
Before melting calculations can be performed, the lithology objects must be assembled in a Mantle object, with their relative mass fractions specified.pyMelt does not limit the number of lithologies that can be assembled in a Mantle object, though in most situations one each of a lherzolite, pyroxenite, and harzburgite lithology are sufficient.The Mantle objects replicate many of the properties of the lithology objects (Figure 1), with methods returning either the mass-weighted average properties, or an array with the value of each lithology.Implicit in our treatment of the lithology objects, and our application of the melting formulation by Phipps Morgan [2001] for adiabatic decompression melting, is an assumption of complete thermal equilibrium but complete chemical disequilibrium between lithologies.

Point melting
The simplest melting calculation that pyMelt can perform is at a single set of conditions; either at a particular  and  using the mantle.F() method, or at a particular  and entropy (expressed by the mantle potential temperature   ) using the mantle.isobaricMelt()method.When calculating the melt fraction for a given  and , pyMelt calls the F() methods for each lithology, returning each in an array.Isentropic melting at a particular  are performed according to the method for isobaric melting in Matthews et al. [2021] in two steps.First the entropy change on cooling the mantle to its solidus is calculated: where   is the   averaged over each lithology.Melting then proceeds by isobaric heating in small temperature increments δ until: where: summing over each lithology , where δ  is the increment of melt fraction corresponding to the increment δ .

Adiabatic decompression melting
Adiabatic decompression melting calculations are performed by the adiabaticMelt method of the Mantle object, requiring only that a value for   is specified.By default the calculation will begin at the solidus and end at 0.01 GPa with a pressure decrement of 0.004 GPa at each decompression step.The calculation proceeds by simultaneously integrating / for each lithology , and / for the melting assemblage, to obtain the melt fractions (  ) of each lithology and the mantle temperature () at each step.The value of / is determined for each melting lithology using Equation 29 of Phipps Morgan [2001] (the mass-weighted average values of   , α, and density ρ are indicated with a bar): The value of / is then obtained from Equation 28 of Phipps Morgan [2001] using the values for one lithology  (arbitrarily, the one with the most negative / in pyMelt): The integration is performed using a fourth-order Runge-Kutta routine.The results of the melting calculation are returned as a meltingColumn object, which records the melt fractions of each lithology, in addition to the aggregate melt fraction and the temperature at each pressure step.
When the calculation is started at a specified pressure, high mantle   may mean the mantle has already exceeded its solidus.In this case an interval of isobaric melting will occur before decompression starts, such that entropy is conserved (Equations 6-8).Decompression melting then proceeds, as described above.
The Phipps Morgan [2001] melting formulation assumes batch melting, whereby the melt is not separated from the solid residue.Therefore, in cases where a lithology  is in thermal equilibrium with a more fusible lithology, the heat extracted by melting of the more fusible lithology can cause lithology  to refreeze, i.e.   / > 0. By default pyMelt will prevent freezing from occurring by setting   / = 0, thereby more closely representing continual melt extraction.This is set as the default behaviour because the chemistry module requires monotonically increasing melt fractions.

Melting hydrous lithologies
pyMelt allows any lithology to be used in decompression melting calculations; however, when hydrous lithologies are modelled by assuming their H 2 O is not extracted it is unlikely to correspond to a realistic melting scenario.Nevertheless, pyMelt can calculate the evolution of melt fraction (including vapour exsolution and the accompanying melt freezing) that such a hydrous lithology would undergo during adiabatic decompression.A demonstration of such calculations is provided in the pyMelt documentation.pyMelt could be extended in the future with new melting methods and geosettings classes (Section 7) to make use of the hydrous lithologies for modelling subduction zone melting.

TRACE ELEMENTS
Following the creation of a meltingColumn object by the adiabaticMelt method, the pyMelt chemistry module can be used to calculate the concentration of each trace element  within the melt  , (Figure 1).The calculation is performed by the meltingColumn.calculateChemistrymethod and requires the concentration of each trace element in each lithology  , , in addition to the parameters required by the chemical model (e.g. the partition coefficients   ).There are four built in chemical models: batch melting, near-fractional melting (instantaneous and accumulated melts), and the INVMEL forward model [McKenzie and O'Nions 1991;White et al. 1992].For the batch and near-fractional melting models the partition coefficient can either be a constant, or a user-defined function of , , and .
Each element (in each lithology) to be included in the calculation is defined as a species object that contains its solid concentration  0 and a composition method for calculating the melt composition as a function of  (and possibly  and ).Defining each element separately permits the incorporation of more complex partitioning behaviour for some elements alongside simpler models for other elements.Generally, users will be unaware of species objects: the meltingColumn.calculateChemistrymethod can assemble them automatically.
The INVMEL model for lherzolite melting incorporates the effects of phase changes and phase exhaustion on the partion-ing of trace elements, in particular the effects of garnet-and clinopyroxene-present melting, but requires many more parameters to be defined.The partition coefficients used by default are those compiled by Gibson and Geist [2010], and other parameters are set to the values used by Ball et al. [2021].
For convenience, the chemistry module has a number of estimates for partition coefficients and mantle trace element concentrations built in (see the documentation for more details).

GEOLOGICAL SETTINGS
In many cases the information provided by the meltingColumn object is sufficient; in other cases a user may be interested in derived properties for a particular geological setting, aggregate melt compositions, for example.The pyMelt geoSetting classes (spreadingCentre and intraPlate) provide this facility, taking a meltingColumn object as an input (Figure 1).
When calculating aggregate properties of the melting region we must consider how each melt in the melting column should be weighted to account for mantle flow.For example, active upwelling in a mantle plume causes more mantle material to pass through the melting region at its base [Maclennan et al. 2001], meaning deeper melts should have a greater weighting in plume models.User-defined weighting functions may be specified for each calculation, but how they are implemented varies between geoSetting classes.An example weighting function built into pyMelt (geosettings.weighting_expdecay)has the form: where λ  and µ are constants, and  max and  min are the maximum and minimum pressures from which melts are formed at the geological setting.

Spreading centres at steady state
When a spreadingCentre object is initialised, the crustal thickness   will be calculated, assuming passive corner-flow mantle upwelling [Plank and Langmuir 1992].To account for the triangular melting region, the total melt fraction is integrated over the melting column, until the pressure exerted by the crust (calculated by stepwise integration with the trapezium rule) is equal to the pressure of the melting step: where  is the acceleration due to gravity on Earth, ρ c is the density of the crust, and  crust is the pressure at the base of the crust.The term  is the optional user-defined weighting function, which takes  = 0 by default.Using the form 1 +  allows separation of the passive upwelling component and the active upwelling component.The (1− φ    ) term in the denominator accounts for compaction, i.e. mantle material will continuously replace the volume lost due to melt extraction [White et al. 1992].The contributions of each lithology to the aggregate crust is calculated similarly: When modelling continental rifts, the pressure exerted by the lithosphere can be imposed and the thickness of igneous crust calculated.
If the meltingColumn object used to generate the spreadingCentre object has chemistry, upon initialisation of the spreadingCentre object, the composition of the homogenised melt is calculated.Similarly to the crustal thickness calculations, homogenisation of chemistry takes into account the triangular melting region, compaction and any additional weighting function.The equation is modified from McKenzie and O'Nions [1991]: and is evaluated using the trapezium rule.If the calculated melt compositions represent instantaneous and not batch melts, each column is first homogenised, with each melt weighted according to its lithology fraction and melt fraction, as above.The weighting function () is not applied in this step.
The crystallisation temperature of melts extracted from the top and base of the melting region can be calculated, using the method described by Matthews et al. [2016].The olivine saturation temperature at the pressure of magma storage is found using the pressure dependence of the olivine saturation surface [39.16K GPa −1 , Putirka 2008].This method is available also in the intraPlate geoSetting class.

Intra-Plate settings at steady state
To initialise an intraPlate geoSetting object, the pressure at the base of the lithosphere must be provided.The calculation results stored in the meltingColumn will be truncated at that pressure.If the relative density of the mantle is provided (Δρ = ρ ambient-mantle − ρ plume-mantle ), the melt flux   is calculated during initialisation using the equation: If the meltingColumn object used to initialise the intraPlate object has chemistry and the melt compositions represent instantaneous melts, they will be homogenised during initialisation according to: If the melt compositions represent accumulated melts and there is a weighting function applied, no result will be returned.

APPLICATIONS
pyMelt has already been used in three published studies, with more studies having been submitted for publication.Matthews et al. [2021] inverted pyMelt to extract estimates of mantle   , pyroxenite fraction, and harzburgite fraction for a number of igneous provinces by matching petrological estimates of crystallisation temperature and geophysically determined magmatic productivity.Gleeson et al. [2021] used pyMelt to estimate the contribution of pyroxenite melts to Galápagos lavas as a function of mantle   and proportion of pyroxenite in the mantle source.Harðardóttir et al. [2022] used pyMelt as the starting point for calculating magma compositions generated by melting a lherzolitic and pyroxenitic mantle, thereby testing models for the petrogenesis of Icelandic basalts.
The implementation of pyMelt in Python enables easy integration with other Python modules, models for running Bayesian inversion calculations for example, as done by Matthews et al. [2021].Since numerical inversion techniques rely on running many forward models, the forward model should be fast.Running an adiabatic decompression melting calculation for a   of 1350 • C and a pure lherzolite mantle takes 455 ms in a Jupyterlab notebook running on a 2019 Mac-Book Pro with 2.6 GHz processor and 16 GB of RAM.A similar calculation, but using a hydrated lherzolite (0.1 wt.% H 2 O modelled as continuous melting) takes much longer (59.6 s on the same machine), due to the need to solve for  numerically at each step.For efficient inversion using a hydrated lithology, further development will be required to increase the computational efficiency.Calculating the trace element composition of the magmas using the default settings takes 7.1 s.
Future development of the library could include methods for calculating melting behaviour at subduction zones, integration with geodynamical models, and other more complex melting scenarios.Work currently underway will provide an expanded set of tools for modelling the behaviour of trace elements during melting of lithologically heterogeneous mantle, and potentially the major element composition of lavas.We welcome contributions to pyMelt from the community and guidelines are available on the GitHub repository.

Figure 1 :
Figure 1: Summary of the structure of pyMelt.Each box represents an instance of a pyMelt object, either created directly by the user (with user-defined variables as specified on the arrows) or returned by a method call.The properties and methods of each class of objects are shown on the right hand side (symbols as defined in the text).
15)which is modified from the equation for volume flux through a deformable conduit[Turcotte and Schubert 2002]. (max) is the melt fraction of lithology  at the top of the conduit,  is the conduit radius (default: 100 km), and µ is the viscosity of the plume (default: 10 19 Pa s), with the default values taken fromShorttle et al. [2014].