Tutorial 1 - Getting Started¶

Simon Matthews (University of Iceland) and Kevin Wong (University of Leeds)

Introduction¶

pyMelt is a python package for calculating the melting behaviour of mantle comprising multiple lithologies. The module implements the melting equations developed by Phipps Morgan (2001) to calculate the melting behaviour of mantle comprising any chosen lithology.

Currently supported calculations:

Adiabatic decompression melting
Isobaric melting

Parameters that can be calculated:

The geotherm for decompressing mantle
Melt fractions for each lithology
Crustal thickness for passive-upwelling at a mid-ocean ridge
Crystallisation temperatures (following the method in Matthews et al., 2016)

Installing pyMelt¶

See the github readme or the readthedocs documentation for installation instructions.

once pyMelt is installed, it can be imported:

import pyMelt as m

adding as m to the end of the instruction allows us to use the shorthand m when accessing the module.

We will also use the numpy, pandas, and matplotlib.pyplot libraries in this tutorial:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

Lithology objects¶

pyMelt offers a number of different lithologies that can be thermodynamically modelled either separately or in combination. pyMelt includes the new parameterisations for KLB-1, KG1, and silica-saturated pyroxenite of Matthews et al., 2021. The lithologies included in the module are:

m.lithologies.matthews.klb1: KLB-1 lherzolite (Matthews et al., 2021)
m.lithologies.matthews.kg1: KG1 silica-undersaturated pyroxenite (Matthews et al., 2021)
m.lithologies.matthews.eclogite: silica-saturated pyroxenite (Matthews et al., 2021)
m.lithologies.shorttle.kg1: KG1 silica-undersaturated pyroxenite (Shorttle et al. 2014)
m.lithologies.katz.lherzolite: lherzolite (Katz et al., 2003)
m.lithologies.pertermann.g2: G2 pyroxenite (Pertermann & Hirschmann, 2002)
m.lithologies.shorttle.harzburgite: non-melting harzburgite

Each lithology is treated as a python object, and an instance can be assigned to a variable:

lz = m.lithologies.matthews.klb1()
px = m.lithologies.matthews.kg1()
hz = m.lithologies.shorttle.harzburgite()

Each lithology object contains methods describing its thermodynamic properties. Most of these methods are hidden and they vary from model to model, depending on how it is formulated. However, every lithology has the following methods:

TSolidus(P): temperature of the lithology solidus (°C) at a given pressure in GPa
TLiquidus(P): temperature of the lithology liquidus (°C) at a given pressure in GPa
F(P, T): melt fraction of the lithology at a given pressure (GPa) and temperature (°C)
dTdF(P, T): dT/dF of the lithology at a constant pressure
dTdP(P, T): dT/dP of the lithology at a constant melt fraction These can be called to get the properties of the lithology at particular temperatures and pressures, for example:

lz_solidus = lz.TSolidus(2.0)
lz_liquidus = lz.TLiquidus(2.0)
lz_F = lz.F(2.0, 1500.0)
print("At 2 GPa the solidus temperature is {:.1f}˚C \n and the liquidus temperature is: {:.1f}˚C. \nAt 1500˚C and 2 GPa the melt fraction is {:.2f}".format(lz_solidus, lz_liquidus, lz_F))

At 2 GPa the solidus temperature is 1397.6˚C 
 and the liquidus temperature is: 1911.3˚C. 
At 1500˚C and 2 GPa the melt fraction is 0.25

Throughout pyMelt the units of temperature are always ˚C, and the units of pressure are always GPa.

We could go one step further and make a plot of the solidus and liquidus temperature for one of the lithologies:

p = np.linspace(0.0,3.0,31)

tliq = np.zeros(np.shape(p))
tsol = np.zeros(np.shape(p))

for i in range(len(p)):
    tliq[i] = lz.TLiquidus(p[i])
    tsol[i] = lz.TSolidus(p[i])

f,a = plt.subplots(dpi=100)
a.plot(p, tliq, label='Liquidus')
a.plot(p, tsol, label='Solidus')
a.set_xlabel('Pressure (GPa)')
a.set_ylabel('Temperature (˚C)')
a.legend()
plt.show()

To see a list of the methods available for any python object, the following can be used within a Jupyter notebook to get access to its documentation.

lz?

Type:        klb1
String form: <pyMelt.lithologies.matthews.klb1 object at 0x7f9973cfe310>
File:        ~/opt/anaconda3/lib/python3.7/site-packages/pyMelt/lithologies/matthews.py
Docstring:  
Implementation of the KLB1 melting model from Matthews et al. (2021).

To use the same format of parameterisation for another lithology, the parameter values
may be changed. They are provided as a dictionary during initialisation of the class.

- Mcpx:   Mass fraction of cpx in the source. Controls the transition to low-productivity
  harzburgite-type melting.
- A1:     Parameter used to define solidus.
- A2:     Parameter used to define solidus.
- A3:     Parameter used to define solidus.
- A4:               Parameter used to define solidus.
- B1:     Parameter used to define liquidus.
- B2:     Parameter used to define liquidus.
- B3:     Parameter used to define liquidus.
- B4:     Parameter used to define liquidus.
- C:        Parameter used to define lherzolite-liquidus.
- beta1:  Parameter used to calculate melt fraction during cpx-present melting.
- beta2:  Parameter used to calculate melt fraction during cpx-absent melting.
- r1:     Parameter used to define cpx reaction coefficient.
- r2:     Parameter used to define cpx reaction coefficient.

The thermal expansivities, the heat capacity, the densities, and the entropy of fusion may
also be changed during class initialisation.

Parameters
----------
CP :         float, default: pyMelt.default_properties['CP']
    The heat capacity (J K-1 kg-1)
alphas :     float, default: pyMelt.default_properties['alphas']
    The thermal expansivity of the solid (1e-6 K-1)
alphaf :     float, default: pyMelt.default_properties['alphaf']
    The thermal expansivity of the melt (1e-6 K-1)
rhos :       float, default: pyMelt.default_properties['rhos']
    The density of the solid (kg m-3)
rhof :       float, default: pyMelt.default_properties['rhof']
    The density of the melt (kg m-3)
DeltaS :     float, default: pyMelt.default_properties['DeltaS']
    The entropy of fusion J K-1 kg-1
parameters : dict, default: parameters from Matthews et al. (2021)

Mantle objects¶

A mantle object is constructed from one or multiple lithologies in specified proportions, and comprises three arguments.

lithologies: a list of the defined lithology objects to be considered in the melting calculation.
proportions: a list of floats (of equivalent length to the list of Argument 1) comprising the relative proportions of the lithologies listed in Argument 1. The floats do not have to be normalised.
names: a list of strings (of equivalent length to the other lists) comprising the names by which the lithologies will be labelled. These strings will be used in data outputs. This can be ommitted, but the results will be unlabelled.

As a demonstration, we can define a three-component mantle. Note that the code is not limited to three lithologies, and can (in principle) have any number of lithologies:

mantle = m.mantle([lz, px, hz], [6, 2, 2], ['Lz', 'Px', 'Hz'])

Here are some of the methods that can be called from the Mantle class:

The bulk thermodynamic properties of the solid mantle can be called:

mantle.bulkProperties()

{'alpha': 40.0, 'CP': 1000.0, 'rho': 3.3}

Does this change once the mantle is molten?

mantle.bulkProperties(T=1400.0, P=2.0)

{'alpha': 42.67799443726366, 'CP': 1000.0, 'rho': 3.261742936610519}

To find out when one of the lithologies will start melting, we can call either the solid_intersection_isobaric method if we are interested in a particular pressure, or we can find out at what pressure the mantle will beginning melting during decompression for a given $T_p$ (potential temperature):

mantle.solidusIntersectionIsobaric(2.0)

array([1397.57332987, 1292.55257242,           inf])

mantle.solidusIntersection(1450.0)

array([3.16884678, 4.45625687,        nan])

Notice that three values were returned, one for each lithology. Since the third lithology is the non-melting harzburgite from Shorttle et al. (2014), it doesn't have a solidus intersection. If we just wanted to know at what temperature the mantle would start melting, but we didn't care which lithology, we could write:

np.nanmin(mantle.solidusIntersectionIsobaric(2.0))

1292.5525724175186

We often talk about mantle temperatures in terms of $T_p$, but sometimes we want to know what temperature this equates to in the solid mantle. The adiabat method will do this. To see the required inputs we can use:

mantle.adiabat?

Signature: mantle.adiabat(P, Tp)
Docstring:
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
-------
float or numpy.array
    Temperature of the mantle at the given pressure and Tp.
File:      ~/opt/anaconda3/lib/python3.7/site-packages/pyMelt/mantle_class.py
Type:      method

So we must specify the pressure and temperature. Note that the documentation specifies this will be the temperature of solid mantle. The adiabatic path taken by semi-molten mantle must be calculated by running a full decompression melting calculation (see below).

mantle.adiabat(10.0, 1300.0)

1502.7229412613215

Stepping back into the world of partially-molten mantle, perhaps we want to know the melt fractions of each lithology at a given temperature and pressure. Just like the pure lithology objects, the F method will perform this calculation:

mantle.F(2.0, 1400.0)

array([0.00185681, 0.47264287, 0.        ])

Adiabatic decompression melting¶

To calculate the consequences of adiabatic decompression melting for this Mantle object, the method AdiabaticMelt can be called which will return a new MeltingColumn object:

column = mantle.adiabaticMelt(1400.0)

This performed the calculation for a $T_p$ of 1400˚C. The way the calculation is performed can be modified, look up the documentation for the other (optional) arguments.

What is going on behind the scenes? This method performs a simultaneous integration of $\frac{dF}{dP}$ and $\frac{dT}{dP}$ to obtain the thermal gradient through the melting region. The melt fraction $F$ of each lithology is then calculated at the same time. Integration is performed using a fourth-order Runge-Kutta algorithm.

pyMelt provides a built in method plot to quickly visualise the results of the calculation.

f,a = column.plot()

Often we want to use these results further, however. We can access the results, like temperature:

column.T

0      1476.497623
1      1476.412794
2      1476.299052
3      1476.180087
4      1476.057254
          ...     
915    1218.612712
916    1218.231565
917    1217.850246
918    1217.468754
919    1217.087089
Name: T, Length: 920, dtype: float64

This returned a pandas.Series object, with 100 values. To keep the notebook tidy, the printout is automatically truncated. But what pressures do these correspond to?

column.P

0      3.688272
1      3.684272
2      3.680272
3      3.676272
4      3.672272
         ...   
915    0.028272
916    0.024272
917    0.020272
918    0.016272
919    0.012272
Name: P, Length: 920, dtype: float64

What about melt fraction?

column.F

0      0.000000
1      0.000038
2      0.000093
3      0.000159
4      0.000231
         ...   
915    0.307799
916    0.308447
917    0.309095
918    0.309744
919    0.310393
Length: 920, dtype: float64

This returns the total melt fraction, but what about the melt fractions for each lithology?

column.lithologies['Lz'].F

0      0.000000
1      0.000000
2      0.000000
3      0.000000
4      0.000000
         ...   
915    0.213295
916    0.213987
917    0.214679
918    0.215372
919    0.216066
Name: F, Length: 920, dtype: float64

column.lithologies['Px'].F

0      0.000000
1      0.000192
2      0.000467
3      0.000793
4      0.001155
         ...   
915    0.899111
916    0.900273
917    0.901436
918    0.902601
919    0.903767
Name: F, Length: 920, dtype: float64

Here we can see the pyroxenite is almost completely molten, but the lherzolite has only achieved 21.5% melting. The pyroxenite has a sufficiently low abundance that it only contributes a small amount to the aggregate melt fraction.

To save this for further work outside of python, we could construct a pandas.DataFrame.

results = pd.DataFrame()

results['T'] = column.T
results['P'] = column.P
results['F_total'] = column.F
results['F_lz'] = column.lithologies['Lz'].F
results['F_px'] = column.lithologies['Px'].F
results

Calling results.to_csv('my_results.csv') will save this as a csv file in your working directory.

Calculating Crustal Thickness at a MOR¶

Crustal thickness can be calculated assuming passive decompression melting in a triangular spreading centre melting region similar to that of a mid-ocean ridge. As part of the pyMelt.geosettings module, there is a spreadingCentre object. When a spreadingCentre object is created from the meltingColumn object, the melt fractions are integrated to obtain the crustal thickness, and more shallow melts are removed.

morb = m.geosettings.spreadingCentre(column)

The new object we created contains the attributes:

tc: integrated crustal thickness at the point where the pressure it exerts is equal to the calculation pressure.
P_base_of_crust: pressure at the base of the crust, at the point where the pressure the generated crust exerts is equal to the calculation pressure.
lithology_contributions integrated proportion of generated crust derived from each lithology where the pressure the generated crust exerts is equal to the calculation pressure. As well as many of the same attributes as the column.

To see the crustal thickness:

morb.tc

11.804219130548248

If we want to know how much each lithology contributes to the crust, we can call the tc_lithology_contributions attribute:

morb.lithology_contributions

{'Lz': 0.18585370789092798, 'Px': 0.814146292109072, 'Hz': 0.0}

In this case the pyroxenite dominates the accumulated melt.

Melt liquidus temperature¶

pyMelt can be used to estimate the liquidus of mantle-derived melts (crystallisation temperature). This is achieved using the MeltCrystallisationT() method. Triangular integration must have been performed beforehand to achieve a liquidus temperature, else an error will be returned.

morb.meltCrystallisationT()

(1251.735708082581, 1347.0844657265302)

This function returns two crystallisation temperature estimates, the first for the melts at the top of the melting column (the shallowest melts), and the second for the melts at the bottom of the melting column (the deepest melts). Both cases assume thermal (and chemical) isolation of the melts during ascent.

See Matthews et al. (2016) or Matthews et al. (2021) for more information about this calculation. To recreate the results from Matthews et al. (2021), the lower crystallisation temperature estimate should be used. Check the documentation to see how the parameters of the calculation may be changed.

Calculating Melt Flux at an Ocean Island¶

Another class in the pyMelt.geosetting module is intraPlate, a very simplistic implementation of a plume upwelling beneath lithosphere. This requires a few more arguments than when creating a spreadingCentre, given by the documentation:

m.geosettings.intraPlate?

Init signature:
m.geosettings.intraPlate(
    MeltingColumn,
    P_lithosphere,
    relative_density=None,
    viscosity=1e+19,
    radius=100000.0,
    weightingFunction=None,
    **kwargs,
)
Docstring:     
Implementation of an intra-plate volcanic province, representing mantle upwelling beneath
lithosphere. The melt flux is calculated assuming flow in a deformable plume conduit (Turcotte
and Schubert, 2002). At present a constant rate of decompression throughout the conduit is
assumed, likely leading to inaccuracies in the estimated melt chemistry.

Parameters
----------
MeltingColumn : pyMelt.meltingcolumn_classes.MeltingColumn
    The melting column from which to construct the geological setting.
P_lithosphere : float, default: 0.0
    The pressure at the base of the lithosphere in a continental rift. If this includes the
    igneous crust, set `extract_melt` to True. Defaults to 0.0, the case of a a mid-ocean
    ridge.
relative_density : float or None, default: None
    The value of (ambient-density - plume-density) in kg m-3.
viscosity : float, default: 1e19
    The viscosity of the mantle plume in Pa s. Default value is 1e19 Pa s, after Shorttle et
    al. (2014).
radius : float, default: 1e5
    The plume radius in m. Default is 1e5 m (or 100 km), after Shorttle et al. (2014).
weighting_function : function or None
    A function of pressure allowing non-uniform weighting of melts throughout the melting
    region. Useful for simulating active upwelling, for example.

Attributes
----------
MeltingColumn : pyMelt.meltingcolumn_classes.MeltingColumn
    The melting column from which the geological setting was constructed
lithologies : dict
    Dictionary containing the DataFrames of the states of each lithology during melting.
mantle : pyMelt.Mantle
    The mantle object from which the melting column was calculated.
P_lithosphere : float
    The pressure at the base of the lithosphere.
melt_flux : float or None
    If the melt flux has been calculated, it will be stored here, in m3 s-1.
lithology_contributions : pandas.Series
    The relative contributions of each lithology to the pooled melt.
chemistry : pandas.Series
    The homogenised melt composition
File:           ~/opt/anaconda3/lib/python3.7/site-packages/pyMelt/geosettings.py
Type:           type
Subclasses:

If you don't wish to calculate a melt flux, the relative_density argument may be ignored. You might want to do this if you are calculating melt chemistry (see tutorial 3), or if you want an easy way of cutting off the top of the melting column. The relative density of the plume compared with ambient mantle determines the upwelling velocity, and therefore the melt production over time. pyMelt does not directly estimate lithology density. The approach used by Shorttle et al. (2014) and Matthews et al. (2021) was to calculate lithology densities as a a function of temperature (at fixed pressure for simplicity) using THERMOCALC. A lookup table was then used to extract an appropriate $\Delta \rho$. To illustrate the calculation we will use an abitrary value of $\Delta \rho = 0.2$ kg m$^{-3}$:

oib = m.geosettings.intraPlate(column, P_lithosphere=1.2, relative_density=0.2)

To print the calculated melt flux (in m$^3$ s$^{-1}$):

oib.melt_flux

1.1609743227613747

	T	P	F_total	F_lz	F_px
0	1476.497623	3.688272	0.000000	0.000000	0.000000
1	1476.412794	3.684272	0.000038	0.000000	0.000192
2	1476.299052	3.680272	0.000093	0.000000	0.000467
3	1476.180087	3.676272	0.000159	0.000000	0.000793
4	1476.057254	3.672272	0.000231	0.000000	0.001155
...	...	...	...	...	...
915	1218.612712	0.028272	0.307799	0.213295	0.899111
916	1218.231565	0.024272	0.308447	0.213987	0.900273
917	1217.850246	0.020272	0.309095	0.214679	0.901436
918	1217.468754	0.016272	0.309744	0.215372	0.902601
919	1217.087089	0.012272	0.310393	0.216066	0.903767