Tutorial 2 - Hydrous Melting

Simon Matthews (University of Iceland), Matthew Gleeson (University of Cardiff), Kevin Wong (University of Leeds)


The lithology objects defined in pyMelt are all anhydrous; however, pyMelt can extend any of these models (and other user defined models) to incorporate the effects of water. We use and extend the formulation set out by Katz et al. (2003), where the solidus temperature is reduced by:

\Delta T(X_{H2O}) = K X^\gamma_{H2O}, \qquad 0 < \gamma < 1

(Eqn. 16 from Katz et al., 2003). X_{H2O} is the melt H_2O content in wt%, and K and \gamma are fitted parameters. The value of X_{H2O} is set either by batch melting (as done by Katz et al., 2003):

X_{H2O} = \frac{X_{H2O}^{bulk}}{D_{H2O} + F (1 - D_{H2O})}

or by continuous melting (i.e., near-fractional melting):

X_{H2O} = \frac{X_{H2O}^{bulk}}{(1-\phi)D + \phi} (1-F)^\frac{(1-\phi)(1-D)}{(1-\phi)D+\phi}

in both equations D is the partition coefficient of H_2O between solid and melt. Additionally, the continous melting equation also requires a value to be set for the porosity \phi, which is the fraction of melt retained and not fractionated.

Note that when using the continuous melting equation it is necessarily not entirely self consistent with the existing melting formulation. The lithology calibrations and the melting equations assume batch melting. However, in most practical applications we make the assumption that this is a close approximation to the melting behaviour of the same lithology during fractional melting. This is clearly not true when water is included- extraction of water will reduce the effect of incorporating water. It is for this reason that modelling water in this manner is still useful.

Katz et al. (2003) limit the amount of water that can enter the melt by imposing a pressure dependent saturation limit. We implement the same equation in pyMelt.

The final modification to the melting equations concerns the calculation of the melt fraction. To calculate \Delta T(X_{H2O}) the value of X_{H2O} (the water content in the melt) must be known. To calculate X_{H2O} the melt fraction is required. There is, therefore, no simple expression to calculate F given all the known parameters. Instead, F must be found by picking a value, calculating the X_{H2O} and F this value imply, and iterating until the picked value of F matches the value calculated.

First, pyMelt must be imported, and other useful packages, must be imported:

import pyMelt
import numpy as np
import matplotlib.pyplot as plt

The hydrous lithology object

The first step is to create a new hydrous lithology object from an anhydrous lithology object. This new object inherits all of the melting behaviour, but with the addition of the above equations (and a suitably adjusted solidus and melt-fraction). We must also specify its H_2O content in wt%. In this case we will create a lithology with 0.1 wt% H_2O, and one with 0.2 wt% H_2O.

To create a hydrated lithology:

# Create an anhydrous lithology object:
lz = pyMelt.lithologies.matthews.klb1()
# Create a hydrated lithology from it:
hlz = pyMelt.hydrousLithology(lz, H2O=0.1)
hhlz = pyMelt.hydrousLithology(lz, H2O=0.2)

Solidus position

We can calculate the new solidii and compare with the anhydrous solidus:

p = np.linspace(0,8.0,100)
t = np.zeros([len(p),4])
for i in range(len(p)):
    t[i,0] = lz.TSolidus(p[i])
    t[i,1] = hlz.TSolidus(p[i])
    t[i,2] = hhlz.TSolidus(p[i])
    t[i,3] = lz.TLiquidus(p[i])

f,a = plt.subplots(dpi=100)
a.plot(t[:,0], p, label='Anhydrous')
a.plot(t[:,1], p, label='0.1 wt%')
a.plot(t[:,2], p, label='0.2 wt%')
a.plot(t[:,3], p, label='Liquidus')
a.set_xlabel('Temperature ($^\circ$C)')
a.set_ylabel('Pressure (GPa)')


The figure shows that the solidus is deplaced to lower temperature by the addition of H_2O, as expected. The liquidus does not change (try changing the calculation above to calculate the liquidus curve from one of the hydrated lithologies). The uptick in solidus temperature at low temperature is due to the limit of H2O saturation in the melt- at lower pressures a vapour phase will become stable in preference to forming a melt.

Melting behaviour (isobaric)

We can also calculate the melt fraction at fixed pressure, in a similar way:

t = np.linspace(900, 2000, 100)
p = 1.5
mf = np.zeros([len(t),4])
for i in range(len(t)):
    mf[i,0] = lz.F(p, t[i])
    mf[i,1] = hlz.F(p, t[i])
    mf[i,2] = hhlz.F(p, t[i])

f,a = plt.subplots(dpi=100)
a.plot(t, mf[:,0], label='Anhydrous')
a.plot(t, mf[:,1], label='0.1 wt%')
a.plot(t, mf[:,2], label='0.2 wt%')
a.set_xlabel('Temperature ($^\circ$C)')
a.set_ylabel('Melt Fraction')


When we defined the lithology objects we didn’t specify which of the two equations for X_{H2O} should be used. In this case pyMelt will default to the batch melting equation. We can see how continuous-melting would compare by defining new lithology objects, with different choices for the parameter \phi (the porosity):

hlz005 = pyMelt.hydrousLithology(lz, 0.1, continuous=True)
hlz001 = pyMelt.hydrousLithology(lz, 0.1, continuous=True, phi=0.001)
hlz01 = pyMelt.hydrousLithology(lz, 0.1, continuous=True, phi=0.01)

And now we can compare the 0.1 wt% melt fraction for both batch and continuous melting regimes:

t = np.linspace(1050, 1400, 100)
p = 1.5
mf = np.zeros([len(t),5])
for i in range(len(t)):
    mf[i,0] = lz.F(p, t[i])
    mf[i,1] = hlz.F(p, t[i])
    mf[i,2] = hlz005.F(p, t[i])
    mf[i,3] = hlz001.F(p, t[i])
    mf[i,4] = hlz01.F(p, t[i])

f,a = plt.subplots(dpi=100)
a.plot(t, mf[:,0], label='Anhydrous', c='0.5')
a.plot(t, mf[:,1], label='0.1 wt% Batch')
a.plot(t, mf[:,2], label='0.1 wt% $\phi$ = 0.5%')
a.plot(t, mf[:,3], label='0.1 wt% $\phi$ = 0.1%')
a.plot(t, mf[:,4], label='0.1 wt% $\phi$ = 1.0%')
a.set_xlabel('Temperature ($^\circ$C)')
a.set_ylabel('Melt Fraction')


Batch vs Continuous makes a very large difference, with the choice of \phi having a much more minor effect.


Just as for the lherzolitic lithology used above, pyMelt offers the opportunity to perform hydrous melting with the pyroxenitic lithologies. Caution must be observed; the formulation by Katz et al. (2003) calibrates the parameters based on hydrous lherzolite melting experiments, adjustment may be required before they are applied to hydrous pyroxenite melting.

Adiabatic Melting

To calculate the melting behaviour of adiabatically decompressing mantle the hydrous lithologies can be treated identically to regular pyMelt lithologies. This means they can form a component of a heterogeneous mantle. Here we will instead consider a pure lithology mantle for the purposes of illustrating the behaviour of the hydrous lithology.

We will use the hydrated lithologies created above to create pyMelt.Mantle objects:

m = pyMelt.mantle([lz], [1.0])
hm = pyMelt.mantle([hlz], [1.0])
hhm = pyMelt.mantle([hhlz], [1.0])
hm005 = pyMelt.mantle([hlz005], [1.0])

hhlz005 = pyMelt.hydrousLithology(lz, 0.2, continuous=True)
hhm005 = pyMelt.mantle([hhlz005], [1.0])

Now run the adiabtic melting calculation for each of these mantle objects (it will take a little while). This time Pstart is set so that melting doesn’t start at a ridiculously high pressure.

c = m.adiabaticMelt(1350.0, Pstart=8.0)
hc = hm.adiabaticMelt(1350.0, Pstart=8.0)
hhc = hhm.adiabaticMelt(1350.0, Pstart=8.0)
hc005 = hm005.adiabaticMelt(1350.0, Pstart=8.0)
hhc005 = hhm005.adiabaticMelt(1350.0, Pstart=8.0)
/Users/simonmatthews/opt/anaconda3/lib/python3.7/site-packages/pyMelt/mantle_class.py:428: UserWarning: Super solidus start
  warn("Super solidus start")

Note that one of the hydrous lithologies reports a super-solidus start. To avoid this the pressure could be increased; however, it is unlikely that the melt fraction is large at the start of the calculation.

f,a = plt.subplots(1, 2, dpi=100, sharey='row')

a[0].plot(c.T, c.P, label='Anhydrous', c='0.7', lw=5)
a[0].plot(hc.T, hc.P, label='0.1 wt%, batch', c='C1')
a[0].plot(hhc.T, hhc.P, label='0.2 wt%, batch', c='C3')
a[0].plot(hc005.T, hc005.P, label='0.1 wt%, cont.', c='C1', ls='--')
a[0].plot(hhc005.T, hhc005.P, label='0.2 wt%, cont.', c='C3', ls='--')

a[1].plot(c.F, c.P, label='Anhydrous', c='0.7', lw=5)
a[1].plot(hc.F, hc.P, label='0.1 wt%, batch', c='C1')
a[1].plot(hhc.F, hhc.P, label='0.2 wt%, batch', c='C3')
a[1].plot(hc005.F, hc005.P, label='0.1 wt%, cont.', c='C1', ls='--')
a[1].plot(hhc005.F, hhc005.P, label='0.2 wt%, cont.', c='C3', ls='--')

a[0].set_xlabel('Temperature ($^\circ$C)')
a[0].set_ylabel('Pressure (GPa)')

As expected, all hydrous models have a low-F melting tail. However, the batch and continuous melting behaviour diverges. For batch melting, an increase in melt fraction persists at most pressures, but for continuous melting the effect rapidly subsides as H_2O is removed from the system. We can see this effect in the crustal thicknesses by establishing SpreadingCentre objects:

morb = pyMelt.geosettings.spreadingCentre(c)
print("Anhydrous: {:.1f} km".format(morb.tc))
morb = pyMelt.geosettings.spreadingCentre(hc)
print("Batch 0.1 wt%: {:.1f} km".format(morb.tc))
morb = pyMelt.geosettings.spreadingCentre(hhc)
print("Batch 0.2 wt%: {:.1f} km".format(morb.tc))
morb = pyMelt.geosettings.spreadingCentre(hc005)
print("Cont. 0.1 wt%: {:.1f} km".format(morb.tc))
morb = pyMelt.geosettings.spreadingCentre(hhc005)
print("Cont. 0.2 wt%: {:.1f} km".format(morb.tc))
Anhydrous: 6.3 km
Batch 0.1 wt%: 9.6 km
Batch 0.2 wt%: 12.3 km
Cont. 0.1 wt%: 6.9 km
Cont. 0.2 wt%: 7.5 km

The increase in crustal thickness due to the presence of H_2O is much more modest for the continuous melting model.

[ ]: