Tutorial 3 - Trace Elements

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

Introduction

pyMelt contains a chemistry module for calculating the trace element composition of melts. There are a number of inbuilt chemical methods: - Batch melting - Continuous melting, which can be used to calculate either instantaneous or accumulated melts - INVMEL, following McKenzie & O’Nions (1991), incorporates the effect of aluminous phase transitions and phase exhaustion Currently the methods all use constant partition coefficients, but soon using pressure and temperature dependent partition coeffients will be an option.

The module is designed to allow maximum customisation of the calculations, but using the default options is very straightforward. INVMEL is used for most elements by default, with the exception of the most incompatible elements, which are modelled assuming continuous melting with 0.5% porosity. See the documentation for more information about the partition coefficients used.

First we must import the required modules:

[1]:
import pyMelt as m
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from copy import copy

Making MORB from a homogeneous mantle

In the first example we will predict the composition of a spreading centre basalt using the default options.

First we must create a mantle object and perform adiabatic melting, as before:

[2]:
lz = m.lithologies.matthews.klb1()
mantle = m.mantle([lz], [1.0], ['lz'])
column = mantle.adiabaticMelt(1350.0)

We must now add chemistry to these melts. The default options are designed for melting depleted mantle, so we don’t need to set up or pass any additional arguments:

[3]:
column.calculateChemistry()
Lithology composition is set to the depleted mantle of Workman & Hart (2005).

Notice that pyMelt told us it was setting the composition to the depleted mantle. See below for how to adjust this. At this point we can already inspect some of the instantaneous results:

[4]:
f, a = plt.subplots(dpi=150)
a.plot(column.P, column.lithologies['lz']['La'], label='La')
a.plot(column.P, column.lithologies['lz']['Yb'], label='Yb')

a.set_ylabel('Concentration (ppmw)')
a.set_xlabel('Pressure (GPa)')
a.legend()
plt.show()
_images/tutorial3_7_0.png

Notice how Yb behaves more compatibly than La, persisting for much longer in melts. The subtle change in gradient in the Yb line occurs during the spinel to plagioclase transition.

When the spreadingCentre geoSetting class is created pyMelt automatically calculates the crustal thickness based on a triangular melting region with passive upwelling (as shown in Tutorial 1). If the MeltingColumn class contains chemistry, creating the geoSetting class also applies the appropriate homogenisation calculation. For a spreading centre this is homogenisation with a triangular weighting.

[5]:
morb = m.geosettings.spreadingCentre(column)

We can view the raw homogenised trace element concentrations:

[6]:
morb.chemistry
[6]:
Rb       0.396402
Ba       4.461680
Th       0.062251
U        0.025427
Nb       1.183375
Ta       0.076496
La       1.509458
Ce       4.252900
Pb       0.142395
Pr       0.807481
Nd       4.248332
Sr      59.324005
Zr      37.916025
Hf       1.073972
Sm       1.648624
Eu       0.654300
Ti    4320.851395
Gd       2.435299
Tb       0.464154
Dy       3.370168
Ho       0.765611
Y       20.813260
Er       2.181567
Yb       2.182753
Lu       0.334373
dtype: float64

But to really make sense of the results it is convenient to make a spider diagram, for which the GeoSetting classes have a built in method. There are many ways to customise the plot, most importantly the reference values used for normalisation. The default normalisation is the primitive mantle of Palme & O’Neill (2013).

[7]:
f,a = morb.plotSpider()
plt.show()
_images/tutorial3_13_0.png

If we wanted to see just the rare earth elements, we can specify this. We can also normalise to chondrite, and show the range of instantaneous melts:

[8]:
f,a = morb.plotSpider(element_order=['La','Ce','Pr','Nd','Sm','Eu','Gd','Tb','Dy','Ho','Er','Yb','Lu'],
                       normalisation='CI', plot_instantaneous=True)
plt.show()
_images/tutorial3_15_0.png

Finally, we can see what the effect of some crystal fractionation would be on the composition. This time I will normalise to depleted mantle, to remove the Pb anomaly. You can also change the axis labels (to specify the normalisation, for example), as shown here.

[9]:
f,a = morb.plotSpider(normalisation='DM', crystal_fraction={'olv':0.05, 'plg':0.2})
a.set_ylabel("Concentration (normalised to DMM)")
plt.show()
_images/tutorial3_17_0.png

Crystallising olivine and plagioclase has enriched the melt in trace elements, with the exception of Sr, since it is very compatible in plagioclase. Eu also shows a negative anomaly owing also to its relative compatibility in plagioclase. The crystallisation calculation can also be used independently of plotting. In fact the plotting method just calls the GeoSetting.crystallisation_chemistry() itself. See the documentation for GeoSetting.crystallisation_chemistry() for information about customising the crystallisation calculation.

Making OIB from a heterogeneous mantle

In addition to calculating the trace element concentrations in melts from a homogeneous mantle, pyMelt allows the trace element contents of melts from multiple lithologies to be calculated simultaneously. However, pyMelt requires additional input to tell it how to treat each individual lithology. This is partly from a practical stand point- one set of default options for pyroxenite (for example) is unlikely to be relevant to many calculations.

Like before we start by setting up a mantle and melting it by adiabatic decompression. We will consider a three component mantle made up of lherzolite, pyroxenite, and non-melting harzburgite.

[10]:
px = m.lithologies.matthews.kg1()
hz = m.lithologies.shorttle.harzburgite()
mantle = m.mantle([lz, px, hz], [0.6, 0.1, 0.3], ['lz','px','hz'])
oi_col = mantle.adiabaticMelt(1480.0)
/Users/simonmatthews/opt/anaconda3/lib/python3.7/site-packages/pyMelt/mantle_class.py:470: UserWarning: Freezing prevented.
  warn("Freezing prevented.")

Notice that pyMelt warned us that freezing was prevented. This can happen when the pyroxenite melting extracts sufficient heat that the partially-molten lherzolite will begin to freeze. When we are using the model to approximate fractional melt extraction we do not want this to happen, and so by default pyMelt will not allow freezing to occur.

pyMelt has some published mantle compositions built in, for example the subducted oceanic crust composition calculated by Stracke et al. (2003):

[11]:
m.chemistry.stracke03_bsic
[11]:
{'Rb': 0.57,
 'Ba': 6.59,
 'Th': 0.088,
 'U': 0.027,
 'Nb': 1.95,
 'Ta': 0.124,
 'La': 1.68,
 'Ce': 5.89,
 'Pb': 0.09,
 'Nd': 7.45,
 'Sr': 81.0,
 'Zr': 64.0,
 'Hf': 1.78,
 'Sm': 2.69,
 'Eu': 1.04,
 'Ti': 7735.0,
 'Gd': 4.03,
 'Dy': 5.01,
 'Y': 28.5,
 'Er': 3.13,
 'Yb': 2.99,
 'Lu': 0.45}

The KG1 bulk composition used as the pyroxenite endmember is a mixture of oceanic crust and lherzolite peridotite. We can mirror this in the trace elements by calculating a mix of the subducted crust composition and the depleted mantle composition. Note that the dict objects have to be converted to pandas.Series objects before applying mathematical operations, and then converted back before feeding to pyMelt.

[12]:
kg1_tes = (0.333 * pd.Series(m.chemistry.stracke03_bsic) +
           0.666 * pd.Series(m.chemistry.palme13_pm))
kg1_tes = kg1_tes[~np.isnan(kg1_tes)]
kg1_tes = dict(kg1_tes)
kg1_tes
[12]:
{'Ba': 6.75657,
 'Ce': 3.1288014,
 'Dy': 2.1504474,
 'Er': 1.3542444,
 'Eu': 0.45720900000000003,
 'Gd': 1.7319330000000002,
 'Hf': 0.7934724000000001,
 'La': 1.0144512,
 'Lu': 0.19702278,
 'Nb': 1.04562,
 'Nd': 3.373956,
 'Pb': 0.15318,
 'Rb': 0.59274,
 'Sm': 1.1852802,
 'Sr': 41.625,
 'Ta': 0.06993,
 'Th': 0.0858474,
 'Ti': 3418.245,
 'U': 0.024242400000000004,
 'Y': 12.24108,
 'Yb': 1.3136184000000002,
 'Zr': 28.1718}

We now have to start setting up the chemistry calculation. We can use the default models again (which for most elements is invmel), but the parameters controlling invmels behaviour will be different for pyroxenite.

To see what parameters are required, we can look at the invmelSpecies object documentation:

[13]:
m.chemistry.invmelSpecies?
Init signature:
m.chemistry.invmelSpecies(
    name,
    c0,
    olv_D,
    cpx_D,
    opx_D,
    spn_D,
    grt_D,
    plg_D,
    mineralProportions=             olv    opx    cpx    grt    spn    plg
grt_field  0.598  0.221  0.076  0.115  0.000  0.000
spn_field  0.578  0.270  0.119  0.000  0.033  0.000
plg_field  0.636  0.263  0.012  0.000  0.000  0.089,
    density=3.3,
    cpxExhaustion=0.18,
    garnetInCoeffs=[666.7, 400.0],
    spinelOutCoeffs=[666.7, 533.0],
    plagioclaseInInterval=[25.0, 35.0],
    **kwargs,
)
Docstring:
Implementation of the forward trace element model used by invmel (McKenzie & O'Nions, 1991).

The default values of the mineral-melt partition coefficients are taken from REFERENCE.

Parameters
----------
name :  str
    The name of the species.
c0 :    float
    The concentration of the species in the solid. Must be the same units as required in the
    output.
olv_D : float
    The partition coefficient for olivine-melt.
cpx_D : float
    The partition coefficient for clinopyroxene-melt.
opx_D : float
    The partition coefficient for orthopyroxene-melt.
spn_D : float
    The partition coefficient for spinel-melt.
grt_D : float
    The partition coefficient for garnet-melt.
plg_D : float
    The partition coefficient for plagioclase-melt.
MineralProportions : pandas.DataFrame, default: pyMelt.chemistry.mo91_MineralProportions
    A dataframe with the proportions of each mineral phase (columns) in the garnet-, spinel-,
    and plagioclase-field for the lithology. See built in defaults for formatting of keys
    and columns: mo91_MineralProportions, klb1_MineralProportions, kg1_MineralProportions.
density : float, default: 3.3
    The density of the mantle (g cm-3)
cpxExhaustion : int, default: 0.18
    The melt fraction at which cpx (and grt/plg/spn) are exhausted.
garnetInCoeffs : list of floats, default: [666.7, 400.0]
    The coefficients to use in the garnet in surface, T = [0] * P - [1]. Default values are
    appropriate for lherzolite.
spinelOutCoeffs : list of floats, default: [666.7, 533.0]
    The coefficients to use in the spinel out surface, T = [0] * P - [1]. Default values are
    appropriate for lherzolite.
plagioclaseInInterval : list of floats, default: [25.0, 35.0]
    The plagioclase in interval in km. Default values are appropriate for lherzolite.

Attributes
----------
name : str
    Name of the species
c0 : float
    Concentration of the species in the solid.
D : numpy.Array
    The partition coefficients in the order: olivine, cpx, opx, spinel, garnet, plag.
MineralProportions_solid : pandas.DataFrame
    See parameters above.
density : float
    The density of the mantle (g cm-3)
File:           ~/opt/anaconda3/lib/python3.7/site-packages/pyMelt/chemistry.py
Type:           type
Subclasses:

Cpx (and grt) persists to much higher melt fractions during pyroxenite melting, though the melt fraction at cpx-exhaustion varies somewhat with pressure. Assuming exhaustion at 60% melting errs towards the low pressure results (it persists for longer at higher pressures)

[14]:
cpxExhaustion = 0.6

For simplicity we will translate the spinel-out reaction from lherzolite to lower pressures for the pyroxenite, and use an identical line for garnet-in:

[15]:
spinelOutCoeffs = [666.7, 50]
garnetInCoeffs = [666.7, 50]

We will not change the plagioclase-spinel interval (it is not so different between KG1 and KLB1).

Perhaps most importantly, we must change the proportions of mineral phases present. Fortunately, pyMelt has the mineral proportions for KG1 stored already in the chemistry module:

[16]:
m.chemistry.kg1_MineralProportions
[16]:
olv opx cpx grt spn plg
grt_field 0.181 0.012 0.422 0.385 0.000 0.000
spn_field 0.110 0.178 0.641 0.000 0.071 0.000
plg_field 0.118 0.150 0.655 0.000 0.000 0.067

We must then feed this information to the chemistry module. Notice that the parameters must be provided for each lithology (even the non-melting harzburgite!). I copied the values for lherzolite from the default values printed in the docstring above.

[17]:
oi_col.calculateChemistry(elements = {'lz': m.chemistry.workman05_dmm,
                                       'px': kg1_tes,
                                       'hz': m.chemistry.workman05_dmm}, # Doesn't matter what value
                           cpxExhaustion = {'lz': 0.18,
                                            'px': cpxExhaustion,
                                            'hz': 0.1}, # Doesn't matter what value
                           garnetInCoeffs = {'lz': [666.7, 400.0],
                                             'px': garnetInCoeffs,
                                             'hz': [0,0]},
                           spinelOutCoeffs = {'lz': [666.7, 533.0],
                                              'px': spinelOutCoeffs,
                                              'hz': [0,0]}, # Doesn't matter what value
                           mineralProportions = {'lz': m.chemistry.klb1_MineralProportions,
                                                 'px': m.chemistry.kg1_MineralProportions,
                                                 'hz': m.chemistry.klb1_MineralProportions} # Doesn't matter what value
                          )

As before, we can view the instantaneous melts within the MeltingColumn object. Here we are comparing the two melting lithologies:

[18]:
f, a = plt.subplots(dpi=150)
a.plot(oi_col.P, oi_col.lithologies['lz']['La'], label='La (lz)', ls='--', c='C0')
a.plot(oi_col.P, oi_col.lithologies['lz']['Yb'], label='Yb (lz)', ls='-', c='C0')

a.plot(oi_col.P, oi_col.lithologies['px']['La'], label='La (px)', ls='--', c='C1')
a.plot(oi_col.P, oi_col.lithologies['px']['Yb'], label='Yb (px)', ls='-', c='C1')

a.set_ylabel('Concentration (ppmw)')
a.set_xlabel('Pressure (GPa)')
a.legend()
plt.show()
_images/tutorial3_36_0.png

To see the homogenised melt that would be erupted through thick lithosphere (1.5 GPa at the base in the example here), we can create an IntraPlate class. We will also apply a greater weighting to the deepest melts in the melting region, to simulate the effect of active upwelling fluxing more mantle through the base of the melting region.

[19]:
oib = m.geosettings.intraPlate(oi_col, 2.5, weightingFunction=m.geosettings.weighting_expdecay,
                               weighting_wavelength=0.1)
[20]:
f,a = oib.plotSpider(element_order=['La','Ce','Nd','Sm','Eu','Gd','Dy','Er','Yb','Lu'])
_images/tutorial3_39_0.png

Changing parameters in the calculation

pyMelt allows you to change any of the parameters used in the calculation. You can set one parameter for every element in every lithology, or set each individually by nesting dictionaries (like the elements variable above). To illustrate this, we could change the olivine partition coefficient for La to make it very compatible.

First, we can obtain the original olivine partition coefficients:

[21]:
ol_D = copy(m.chemistry.olv_D)
ol_D
[21]:
{'Rb': 0.0003,
 'Ba': 5e-06,
 'Th': 5e-05,
 'U': 0.00038,
 'Nb': 0.0005,
 'Ta': 0.0005,
 'La': 0.0005,
 'Ce': 0.0005,
 'Pb': 0.003,
 'Pr': 0.0008,
 'Nd': 0.00042,
 'Sr': 4e-05,
 'Zr': 0.0033,
 'Hf': 0.0022,
 'Sm': 0.0011,
 'Eu': 0.0016,
 'Ti': 0.015,
 'Gd': 0.0011,
 'Tb': 0.0015,
 'Dy': 0.0027,
 'Ho': 0.0016,
 'Y': 0.0099,
 'Er': 0.013,
 'Yb': 0.02,
 'Lu': 0.02}

We must copy the variable otherwise changing the La content directly will change the default for every calculation. Now we can selectively change the La partition coefficient:

[22]:
ol_D['La'] = 2.0

To repeat the calculation above, but with the new olivine partition coefficients affecting all the lithologies:

[23]:
oi_col.calculateChemistry(olv_D = ol_D,
                           elements = {'lz': m.chemistry.workman05_dmm,
                                       'px': kg1_tes,
                                       'hz': m.chemistry.workman05_dmm}, # Doesn't matter what value
                           cpxExhaustion = {'lz': 0.18,
                                            'px': cpxExhaustion,
                                            'hz': 0.1}, # Doesn't matter what value
                           garnetInCoeffs = {'lz': [666.7, 400.0],
                                             'px': garnetInCoeffs,
                                             'hz': [0,0]},
                           spinelOutCoeffs = {'lz': [666.7, 533.0],
                                              'px': spinelOutCoeffs,
                                              'hz': [0,0]}, # Doesn't matter what value
                           mineralProportions = {'lz': m.chemistry.klb1_MineralProportions,
                                                 'px': m.chemistry.kg1_MineralProportions,
                                                 'hz': m.chemistry.klb1_MineralProportions} # Doesn't matter what value
                          )
/Users/simonmatthews/opt/anaconda3/lib/python3.7/site-packages/pyMelt/chemistry.py:689: UserWarning: Discretisation is too course to capture the behaviour of La.
  warn("Discretisation is too course to capture the behaviour of " + self.name + ".")
[24]:
oib_new = m.geosettings.intraPlate(oi_col, 2.5, weightingFunction=m.geosettings.weighting_expdecay,
                               weighting_wavelength=0.1)
oib_new.plotSpider(element_order=['La','Ce','Nd','Sm','Eu','Gd','Dy','Er','Yb','Lu'])
[24]:
(<Figure size 900x600 with 1 Axes>,
 <matplotlib.axes._subplots.AxesSubplot at 0x7fb3b6087e50>)
_images/tutorial3_46_1.png

If we instead wanted to make every element in garnet incompatible, we could pass a single float:

[25]:
oi_col.calculateChemistry(grt_D = 1e-3,
                           elements = {'lz': m.chemistry.workman05_dmm,
                                       'px': kg1_tes,
                                       'hz': m.chemistry.workman05_dmm}, # Doesn't matter what value
                           cpxExhaustion = {'lz': 0.18,
                                            'px': cpxExhaustion,
                                            'hz': 0.1}, # Doesn't matter what value
                           garnetInCoeffs = {'lz': [666.7, 400.0],
                                             'px': garnetInCoeffs,
                                             'hz': [0,0]},
                           spinelOutCoeffs = {'lz': [666.7, 533.0],
                                              'px': spinelOutCoeffs,
                                              'hz': [0,0]}, # Doesn't matter what value
                           mineralProportions = {'lz': m.chemistry.klb1_MineralProportions,
                                                 'px': m.chemistry.kg1_MineralProportions,
                                                 'hz': m.chemistry.klb1_MineralProportions} # Doesn't matter what value
                          )
[26]:
oib_new = m.geosettings.intraPlate(oi_col, 2.5, weightingFunction=m.geosettings.weighting_expdecay,
                                   weighting_wavelength=0.1)
oib_new.plotSpider(element_order=['La','Ce','Nd','Sm','Eu','Gd','Dy','Er','Yb','Lu'])
[26]:
(<Figure size 900x600 with 1 Axes>,
 <matplotlib.axes._subplots.AxesSubplot at 0x7fb3b6797810>)
_images/tutorial3_49_1.png

This is a some what useless and unlikely use case, but universally setting parameters could be useful in other contexts. For example, if using the continuous melting model you might want to set the porosity during melting universally.

Using a variable partition coefficient

The species objects that are built into pyMelt can incorporate partition coefficients that vary with pressure and temperature (for example). To demonstrate this we can use the partitioning of Na into cpx, which increases substantially with pressure.

First we have to create a function that takes in pressure. This function approximates the line shown in Figure 9 of Jennings & Holland (2015) for the bulk partition coefficient.

[31]:
def DNa(state):
    return 0.03*np.exp(state['P']/1.8)

We can demonstrate exactly what this looks like:

[32]:
f,a = plt.subplots()

x = np.linspace(0,6,10)
y = np.zeros(np.shape(x))

for i in range(len(x)):
    y[i] = DNa({'P':x[i]})

a.plot(x, y)

a.set_yscale('log')
a.set_ylim(0.01,1)

a.set_xlabel('Pressure (GPa)')
a.set_ylabel('D(cpx)')

plt.show()
_images/tutorial3_54_0.png

This can be passed in a dictionary in place of a numerical value of the partition coefficient, but if we want to calculate just one element it may be faster to create the species object directly:

[33]:
species = m.chemistry.continuousSpecies_instantaneous('Na', 0.13, DNa)

Setting up the problem like this means that the calculations will only be accurate for cpx-present melting. We can see how this affects melts as a function of pressure by iterating through some different mantle temperatures:

[34]:
Tps = [1350, 1400, 1450, 1500, 1550]

lz = m.lithologies.matthews.klb1()
mantle = m.mantle([lz], [1.0], ['lz'])

f,a = plt.subplots(dpi=150)

for i in range(len(Tps)):
    column = mantle.adiabaticMelt(Tps[i])
    column.calculateChemistry(species_objects = {'lz' :[species]})
    a.plot(column.P, column.lithologies['lz']['Na'], label='Tp = {:.0f}$^\circ$C'.format(Tps[i]))

a.legend()
a.set_xlabel('Pressure (GPa)')
a.set_ylabel('Na$_2$O (wt%)')

plt.show()
_images/tutorial3_58_0.png
[ ]:

[ ]:

[ ]: