pyMelt package

Subpackages

Submodules

pyMelt.core module

exception pyMelt.core.ConvergenceError(message)[source]

Bases: Error

Exception raised for errors in converging.

Attributes:

expression – input expression in which the error occurred message – explanation of the error

exception pyMelt.core.Error[source]

Bases: Exception

Base class for exceptions in this module.

exception pyMelt.core.InputError(message)[source]

Bases: Error

Exception raised for errors in the input.

Attributes:

expression – input expression in which the error occurred message – explanation of the error

pyMelt.lithology_class module

pyMelt.mantle_class module

class pyMelt.mantle_class.mantle(lithologies, proportions, names=None)[source]

Bases: object

The mantle class consists of one or more lithology objects, in a particular proportion. The mantle class contains the methods used for doing melting calculations and for calculating the properties of a heterogeneous mantle.

Parameters:

  • lithologies (list of lithology objects) – A list of lithology instances.

  • proportions (list or numpy.array of floats) – The mass ratios of the lithologies, doesn’t need to be normalised.

  • names (list of strings or None, default: None) – The names of the lithologies. If False, default names will be chosen.

Variables:

  • number_lithologies (int) – the number of lithologies in the mantle class

  • CP (list of floats) – the heat capacities of the lithologies

  • alphaf (list of floats) – the thermal expansion coefficients of the melts produced by each lithology.

  • alphas (list of floats) – the thermal expansion coefficients of each lithology

  • rhof (list of floats) – the densities of the melts produced by each lithology

  • rhos (list of floats) – the densities of each lithology

  • DeltaS (list of floats) – the entropy change on melting of each lithology.

F(P, T)[source]

Calculates the melt fraction of each lithology at a given pressure and temperature.

Parameters:

  • P (float) – Pressure in GPa

  • T (float) – Temperature in degC

Returns:

Array containing the melt fraction of each lithology.

Return type:

numpy.Array

adiabat(P, Tp)[source]

Calculates the actual temperature of the solid mantle at a given pressure, given the potential temperature.

Parameters:

  • P (float or numpy.array) – Pressure in GPa.

  • Tp (float or numpy.array) – Potential temperature in degC.

Returns:

Temperature of the mantle at the given pressure and Tp.

Return type:

float or numpy.array

adiabaticGradient(P, T)[source]

Calculates dTdP if melting has gone to completion (or hasn’t started) for the bulk mantle.

Parameters:

  • P (float) – Pressure in GPa.

  • T (float) – Temperature in degC.

Returns:

The adiabatic gradient in C/GPa

Return type:

float

adiabaticMelt(Tp, Pstart=None, Pend=0.01, dP=-0.004, steps=None, ReportSSS=True, adjust_pressure=True, prevent_freezing=True, warn_prevent_freezing=True)[source]

Performs simultaneous integration of dFdP and dTdP to obtain the thermal gradient through the melting region. F of each lithology is then calculated using the P,T path. Integration is performed using a 4th order Runge-Kutta algorithm.

The T-P path is allowed to overstep the solidus on the step prior to the start of melting.

Parameters:

  • Tp (float) – The potential temperature (degC) at which to perform the calculation.

  • Pstart (float or None, default: None) – The pressure (in GPa) at which to begin upwelling. If None, the calculation will start at the solidus.

  • Pend (float, default: 0.0) – The pressure (in GPa) at which to stop upwelling.

  • dP (float, default: -0.004) – The step size in pressure (GPa) to use in the calculation. If the argument steps is set, dP will be ignored.

  • steps (None or int, default: None) – The number of dP increments to split the melting region into. If set to None, the number of steps is determined by dP.

  • ReportSSS (bool, default: True) – Print to the console if the start is above the solidus of one of the lithologies. Either way the code will calculate the melt fraction at this point by conserving entropy. Set to False if you don’t want to be warned.

  • adjust_pressure (bool, default: True) – If True, the pressure range will be adjusted slightly so that one point coincides with the solidus. This should avoid issues with discretisation.

  • prevent_freezing (bool, default: True) – In some melting regions heat extraction by one lithology can cause another’s melts to partially freeze (releasing some heat). If set to True this is prevented from happening. This is most useful when modelling fractional melt extraction.

  • warn_prevent_freezing (bool, default: True) – If set to True, when a melt is prevented from freezing, a warning will be raised.

Returns:

The results are returned in a 1D Melting Column instance, further calculations, e.g. crustal thickness may then be performed on this instance, if desired.

Return type:

pyMelt.meltingColumn

bulkProperties(P=None, T=None)[source]

Calculates the bulk thermodynamic properties of the solid or partially molten mantle.

Parameters:

  • P (float or None, default: None) – The pressure of interest. If None, the properties of the solid mantle will be returned.

  • T (float or None, default: None) – The temperature of interest. If None, the properties of the solid mantle will be returned.

Returns:

The bulk alpha, CP and rho for the mantle at the given P and T, labelled as such.

Return type:

dict

dFdP(P, T, prevent_freezing=False, F_prev=None)[source]

Calculates the value of dFdP for each lithology at the given pressure and temperature, using Eq(26) of Phipps Morgan (2001).

Parameters:

  • P (float) – Pressure in GPa.

  • T (float) – Temperature in degC.

  • prevent_freezing (bool, default: False) – If set to True, any dFdP values > 0 will be set to 0.

  • F_prev (numpy.Array or None, default:None) – If preventing freezing, this is the melt fraction on the previous step to use when checking for freezing.

Returns:

Array of dFdP values for each lithology.

Return type:

numpy.Array

dTdP(P, T, dFdP=None)[source]

Calculates dTdP using Eq(28) of Phipps Morgan (2001). Picks the lithology to use by the one with the largest increase in melt fraction with decompression (though this choice shouldn’t matter).

Parameters:

  • P (float) – Pressure in GPa.

  • T (float) – Temperature in degC.

  • dFdP (numpy.Array or None, default: None) – The value of dFdP at the same T and P. If None, the dFdP method will be called. In most melting calculations the value will already have been calculated, so passing the value will save computational time.

Returns:

The thermal gradient in the melting region at the P and T of interest.

Return type:

float

isobaricMelt(Tstart, P, dT=0.1)[source]

Calculates the amount of melt generated, and the mantle’s temperature, after an interval of melting occuring due to mantle being instantaneously placed above its solidus.

The intention of this function was to handle high Tp cases where the solidus is always exceeded, not to produce physically meaningful results, but to allow the tails of Tp distributions to be approximated reasonably when inverting.

Parameters:

  • Tstart (float) – The temperature (degC) at which to place the solid mantle.

  • P (float) – The pressure at which to perform the calculation.

  • dT (float) – The interval of discretisation of temperature increase from the solidus.

Returns:

Temperature of mantle following the melting step.

Return type:

float

solidusIntersection(Tp)[source]

Finds the pressure at which each lithology’s solidus will be intersected, assuming the mantle follows the solid adiabat up until that point.

Parameters:

Tp (float) – The mantle potential temperature in degC.

Returns:

The pressure of solidus intersection of each lithology.

Return type:

numpy.array

solidusIntersectionIsobaric(P)[source]

Finds the pressure at which each lithology’s solidus will be intersected, assuming the mantle is heated isobarically.

Parameters:

P (loat) – The pressure of interest in GPa

Returns:

The temperature of solidus intersection of each lithology.

Return type:

numpy.Array

pyMelt.meltingcolumn_classes module

Melting Columns

the meltingcolumn_classes module provides the melting column classes. At the moment it consists of a single melting column- a simple 1D melting column.

class pyMelt.meltingcolumn_classes.meltingColumn(calculation_results, mantle, Tp=None)[source]

Bases: object

Class for storing the results of a 1D multi-lithology melting model.

Parameters:

  • calculation_results (pandas.DataFrame) – Dataframe with columns ‘P’ for Pressure in GPa, ‘T’ for Temperature in degrees C, Remaining columns for melt fraction from each lithology.

  • mantle (pyMelt.Mantle) – The mantle object used to generate the melting column.

  • Tp (float) – The potential temperature used to generate the melting column, if applicable.

Variables:

  • calculation_results (pandas.DataFrame) – The stored raw calculation results.

  • mantle (pyMelt.Mantle) – The mantle object used to generate the melting column.

  • Tp (float) – The potential temperature used to generate the melting column, if applicable.

  • P (pandas.Series) – The pressure steps in GPa

  • T (pandas.Series) – The temperature steps in degC.

calculateChemistry(elements=None, species_objects=None, method='default', **kwargs)[source]

Calculate the composition of the melt according to default (or user defined) chemical models. The method may be run with default options if the mantle consists of only one lithology. Otherwise the parameters for each lithology must be specified, or pre-defined species objects must be provided.

If elements is not set and the mantle is made of one lithology, the composition will be set to the depleted mantle of Workman & Hart (2005). The INVMEL forward model is used by default, except for Ba and Rb which are modelled using continuous_instantaneous. The mineral-specific partition coefficients for INVMEL are the constant values compiled by Gibson & Geist (2010). The bulk partition coefficients for Ba and Rb are from Workman & Hart (2005).

Unless passing the species objects directly, the parameters used by the trace element model must be given also. If a constant value should be used for every lithology and every element (e.g., phi for continuous melting) it can be passed as a float. If a value varies for each element, but not each lithology, for example a mineral-specific D, a dict can be passed with each element name as a key. If the parameter varies for each element and lithology, it can be supplied as a nested-dictionary, much like the elements parameter. See Notes for a list of parameters for the default models.

Parameters:

  • elements (dict of dicts of floats (optional)) – A dictionary with each lithology as a key, the values are themselves dictionaries with species names as keys, and their concentrations as the values. The elements must be among the default element list. Ignore if supplying the species objects directly in the species_objects argument.

  • species_objects (dict or None, default: None) – Each lithology name is a key of the dictionary, the values are a list of pyMelt.chemical_classes.species objects. If None, the species objects will be generated from the elements argument.

  • method (string or dict or dict of dicts, default: ‘default’) – One of ‘default’, ‘batch’, ‘continuous_accumulated’, ‘continuous_instantaneous’, ‘invmel’. If using different models for different elements, specify them as a dictionary. This can be nested within another dictionary if you wish to use different combinations for each lithology.

Notes

The ‘batch’ melting routine uses:

  • D, the bulk partition coefficient

The ‘continuous’ melting routine uses:

  • D, the bulk partition coefficient

  • phi, the porosity during melting (as a fraction, not percent).

The ‘invmel’ melting routine uses:

  • olv_D, the olivine-melt partition coefficient

  • cpx_D, the olivine-melt partition coefficient

  • opx_D, the olivine-melt partition coefficient

  • spn_D, the olivine-melt partition coefficient

  • grt_D, the olivine-melt partition coefficient

  • plg_D, the olivine-melt partition coefficient

  • MineralProportions, the mass fraction of each mineral present (prior to melting) in the spinel, plagioclase, and garnet field.

  • cpxExhaustion, the melt fraction at which cpx (and grt/plg/spn) are exhausted.

  • garnetInCoeffs, coefficients controlling the P and T of the garnet in reaction

  • spinelOutCoeffs, coefficients controlling the P and T of the spinel out reaction

  • plagioclaseInInterval, The plagioclase in interval (in km).

plot()[source]

Generates a plot showing the thermal gradient and melt fractions of each lithology.

Returns:

The generated figure and axes objects.

Return type:

(matplotlib.figure, matplotlib.axes)

Module contents

pyMelt

A python package for performing mantle melting calculations.

Developed and maintained by Simon Matthews (simonm@hi.is), Kevin Wong, and Matthew Gleeson.

This module facilitates calculations for melting of heterogeneous mantle, based on empirical parameterisations of the melting functions. The calculations assume complete thermal equilibrium and complete chemical disequilibrium between adjacent mantle lithologies. Various tools are provided for performing decompression melting calculations.

class pyMelt.hydrousLithology(lithology, H2O, D=0.01, K=43.0, gamma=0.75, chi1=12.0, chi2=1.0, l=0.6, continuous=False, phi=0.005, threshold_F=None)[source]

Bases: object

The hydrous lithology class modifies the melting expressions in a given lithology object so that water-present melting can be modelled.

Parameters:

  • lithology (pyMelt.lithology_classes.lithology instance) – The anhydrous lithology for which hydrous melting should be modelled.

  • H2O (float) – The water content of the lithology in wt%.

  • D (float, default: 0.01) – Partition coefficient for water partitioning between melt and the (bulk) residue.

  • K (float, default: 43.0) – Parameter controlling hydrous solidus position, from Katz et al. (2003). See notes.

  • gamma (float, default: 0.75) – Parameter controlling the hydrous solidus position, from Katz et al. (2003). Must take a value between 0 and 1. See notes.

  • chi1 (float, default: 12.0) – Parameter controlling the H2O content of melt at H2O-satuation, in wt% GPa^(-l).

  • chi2 (float, default: 1.0) – Parameter controlling the H2O content of melt at H2O-satuation, in wt% GPa^(-1).

  • l (float, default: 0.6) – Parameter controlling the H2O content of melt at H2O-satuation.

  • continuous (bool, default: False) – Controls whether water extracted to melt is kept or removed from the system. If False, melting will be done assuming batch closed-system melting. If True, the water contents are modelled according to continuous melting, i.e., near-fractional melting, with the porosity controlled by the phi argument.

  • phi (float, default: 0.005) – The porosity to use during continuous melting, if applicable. The default value is 0.5%.

  • threshold_F (float or None, default: None) – If set, when the melt fraction goes above the chosen threshold, the hydrous extension will be skipped. Doing so may help speed up the calculations, but a threshold melt fraction should be chosen with care.

Notes

The parameters controlling the position of the hydrous solidus are defined by eqn. 16 from Katz et al. (2003):

\Delta T(X_{H2O}) = K X^\gamma_{H2O}

Where X_{H2O} is the water content in the melt.

The water content at saturation is controlled by eqn. 17 of Katz et al. (2003):

X^{sat}_{H2O} = \chi_1 P^\lambda + \chi_2 P

F(P, T, **kwargs)[source]

Returns the melt fraction of the hydrated lithology, by iteratively finding the value of F which provides the amount of water to the melt such that the same value of F is returned by the calculation. This is equivalent to Eqn 19 of Katz et al. (2003), but its form will depend on the lithology being used.

The numerical problem is solved using the SciPy method root_scalar from the optimize module. The method used is brentq, so that the problem may be constrained between a melt fraction of 0 and 1.

Parameters:

  • P (float) – Pressure in GPa

  • T (float) – Temperature in degC.

Returns:

The melt fraction.

Return type:

float

H2O_saturation(P)[source]

The H2O content of the melt at H2O saturation. This is Eqn 17 of Katz et al. (2003).

Parameters:

P (float) – Pressure in GPa

Returns:

The H2O content of the melt at H2O-saturation, in wt%.

Return type:

float

TLiquidus(P, **kwargs)[source]

Returns the temperature of the liquidus at any given pressure.

Parameters:

P (float) – Pressure (GPa).

Returns:

Liquidus temperature (degC).

Return type:

float

TSolidus(P, F=0.0)[source]

Returns the temperature of the solidus at any given pressure. Since the solidus is a function of the water content of the bulk, the effective solidus position will shift during melting.

Parameters:

  • P (float) – Pressure (GPa).

  • F (float, default: 0.0) – Melt fraction, between 0 and 1.

Returns:

Solidus temperature (degC).

Return type:

float

dTdF(P, T, **kwargs)[source]

Calculates dT/dF(const. P) for the hydrated lithology numerically.

Parameters:

  • P (float) – Pressure (GPa)

  • T (float) – Temperature (degC)

Returns:

dT/dF(const. P) (K).

Return type:

float

dTdP(P, T, **kwargs)[source]

Calculates dT/dP(const. F) for the hydrated lithology numerically.

The method uses the root_scalar method of the scipy.optimize module to find the curve of T-P for which F=const. The scipy derivative method is then used to numerically differentiate this curve.

Parameters:

  • P (float) – Pressure (GPa).

  • T (float) – Temperature (degC).

Returns:

dT/dP(const. F) (K GPa-1).

Return type:

float

melt_H2O(F, P=None)[source]

Returns the H2O content of the melt at a particular melt fraction, assuming either batch or continuous melting, and H2O saturation controlled by the H2O_saturation method. If H2O is saturated, it will remain in the bulk system.

Parameters:

  • F (float) – The melt fraction. Must be between 0 and 1.

  • P (float or None, default: None) – The pressure in GPa. Used for checking for H2O saturation. If set to None this check will be bypassed.

Returns:

The H2O content of the melt, in wt%.

Return type:

float

class pyMelt.mantle(lithologies, proportions, names=None)[source]

Bases: object

The mantle class consists of one or more lithology objects, in a particular proportion. The mantle class contains the methods used for doing melting calculations and for calculating the properties of a heterogeneous mantle.

Parameters:

  • lithologies (list of lithology objects) – A list of lithology instances.

  • proportions (list or numpy.array of floats) – The mass ratios of the lithologies, doesn’t need to be normalised.

  • names (list of strings or None, default: None) – The names of the lithologies. If False, default names will be chosen.

Variables:

  • number_lithologies (int) – the number of lithologies in the mantle class

  • CP (list of floats) – the heat capacities of the lithologies

  • alphaf (list of floats) – the thermal expansion coefficients of the melts produced by each lithology.

  • alphas (list of floats) – the thermal expansion coefficients of each lithology

  • rhof (list of floats) – the densities of the melts produced by each lithology

  • rhos (list of floats) – the densities of each lithology

  • DeltaS (list of floats) – the entropy change on melting of each lithology.

F(P, T)[source]

Calculates the melt fraction of each lithology at a given pressure and temperature.

Parameters:

  • P (float) – Pressure in GPa

  • T (float) – Temperature in degC

Returns:

Array containing the melt fraction of each lithology.

Return type:

numpy.Array

adiabat(P, Tp)[source]

Calculates the actual temperature of the solid mantle at a given pressure, given the potential temperature.

Parameters:

  • P (float or numpy.array) – Pressure in GPa.

  • Tp (float or numpy.array) – Potential temperature in degC.

Returns:

Temperature of the mantle at the given pressure and Tp.

Return type:

float or numpy.array

adiabaticGradient(P, T)[source]

Calculates dTdP if melting has gone to completion (or hasn’t started) for the bulk mantle.

Parameters:

  • P (float) – Pressure in GPa.

  • T (float) – Temperature in degC.

Returns:

The adiabatic gradient in C/GPa

Return type:

float

adiabaticMelt(Tp, Pstart=None, Pend=0.01, dP=-0.004, steps=None, ReportSSS=True, adjust_pressure=True, prevent_freezing=True, warn_prevent_freezing=True)[source]

Performs simultaneous integration of dFdP and dTdP to obtain the thermal gradient through the melting region. F of each lithology is then calculated using the P,T path. Integration is performed using a 4th order Runge-Kutta algorithm.

The T-P path is allowed to overstep the solidus on the step prior to the start of melting.

Parameters:

  • Tp (float) – The potential temperature (degC) at which to perform the calculation.

  • Pstart (float or None, default: None) – The pressure (in GPa) at which to begin upwelling. If None, the calculation will start at the solidus.

  • Pend (float, default: 0.0) – The pressure (in GPa) at which to stop upwelling.

  • dP (float, default: -0.004) – The step size in pressure (GPa) to use in the calculation. If the argument steps is set, dP will be ignored.

  • steps (None or int, default: None) – The number of dP increments to split the melting region into. If set to None, the number of steps is determined by dP.

  • ReportSSS (bool, default: True) – Print to the console if the start is above the solidus of one of the lithologies. Either way the code will calculate the melt fraction at this point by conserving entropy. Set to False if you don’t want to be warned.

  • adjust_pressure (bool, default: True) – If True, the pressure range will be adjusted slightly so that one point coincides with the solidus. This should avoid issues with discretisation.

  • prevent_freezing (bool, default: True) – In some melting regions heat extraction by one lithology can cause another’s melts to partially freeze (releasing some heat). If set to True this is prevented from happening. This is most useful when modelling fractional melt extraction.

  • warn_prevent_freezing (bool, default: True) – If set to True, when a melt is prevented from freezing, a warning will be raised.

Returns:

The results are returned in a 1D Melting Column instance, further calculations, e.g. crustal thickness may then be performed on this instance, if desired.

Return type:

pyMelt.meltingColumn

bulkProperties(P=None, T=None)[source]

Calculates the bulk thermodynamic properties of the solid or partially molten mantle.

Parameters:

  • P (float or None, default: None) – The pressure of interest. If None, the properties of the solid mantle will be returned.

  • T (float or None, default: None) – The temperature of interest. If None, the properties of the solid mantle will be returned.

Returns:

The bulk alpha, CP and rho for the mantle at the given P and T, labelled as such.

Return type:

dict

dFdP(P, T, prevent_freezing=False, F_prev=None)[source]

Calculates the value of dFdP for each lithology at the given pressure and temperature, using Eq(26) of Phipps Morgan (2001).

Parameters:

  • P (float) – Pressure in GPa.

  • T (float) – Temperature in degC.

  • prevent_freezing (bool, default: False) – If set to True, any dFdP values > 0 will be set to 0.

  • F_prev (numpy.Array or None, default:None) – If preventing freezing, this is the melt fraction on the previous step to use when checking for freezing.

Returns:

Array of dFdP values for each lithology.

Return type:

numpy.Array

dTdP(P, T, dFdP=None)[source]

Calculates dTdP using Eq(28) of Phipps Morgan (2001). Picks the lithology to use by the one with the largest increase in melt fraction with decompression (though this choice shouldn’t matter).

Parameters:

  • P (float) – Pressure in GPa.

  • T (float) – Temperature in degC.

  • dFdP (numpy.Array or None, default: None) – The value of dFdP at the same T and P. If None, the dFdP method will be called. In most melting calculations the value will already have been calculated, so passing the value will save computational time.

Returns:

The thermal gradient in the melting region at the P and T of interest.

Return type:

float

isobaricMelt(Tstart, P, dT=0.1)[source]

Calculates the amount of melt generated, and the mantle’s temperature, after an interval of melting occuring due to mantle being instantaneously placed above its solidus.

The intention of this function was to handle high Tp cases where the solidus is always exceeded, not to produce physically meaningful results, but to allow the tails of Tp distributions to be approximated reasonably when inverting.

Parameters:

  • Tstart (float) – The temperature (degC) at which to place the solid mantle.

  • P (float) – The pressure at which to perform the calculation.

  • dT (float) – The interval of discretisation of temperature increase from the solidus.

Returns:

Temperature of mantle following the melting step.

Return type:

float

solidusIntersection(Tp)[source]

Finds the pressure at which each lithology’s solidus will be intersected, assuming the mantle follows the solid adiabat up until that point.

Parameters:

Tp (float) – The mantle potential temperature in degC.

Returns:

The pressure of solidus intersection of each lithology.

Return type:

numpy.array

solidusIntersectionIsobaric(P)[source]

Finds the pressure at which each lithology’s solidus will be intersected, assuming the mantle is heated isobarically.

Parameters:

P (loat) – The pressure of interest in GPa

Returns:

The temperature of solidus intersection of each lithology.

Return type:

numpy.Array