import numpy as _np
from warnings import warn as _warn
import pandas as _pd
from scipy.optimize import root_scalar as _root_scalar
from pyMelt.meltingcolumn_classes import meltingColumn as _meltingColumn
from pyMelt.core import InputError
[docs]class mantle:
"""
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.
Attributes
----------
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.
"""
def __init__(self, lithologies, proportions, names=None):
self.lithologies = lithologies
if isinstance(proportions, list):
proportions = _np.array(proportions)
self.proportions = proportions / sum(proportions)
self.number_lithologies = _np.shape(self.lithologies)[0]
if names is None:
names = list()
for i in range(self.number_lithologies):
names.append('default ' + str(i))
self.names = names
self.CP = _np.zeros(self.number_lithologies)
self.alphaf = _np.zeros(self.number_lithologies)
self.alphas = _np.zeros(self.number_lithologies)
self.rhof = _np.zeros(self.number_lithologies)
self.rhos = _np.zeros(self.number_lithologies)
self.DeltaS = _np.zeros(self.number_lithologies)
for i in range(self.number_lithologies):
self.CP[i] = self.lithologies[i].CP
self.alphaf[i] = self.lithologies[i].alphaf
self.alphas[i] = self.lithologies[i].alphas
self.rhof[i] = self.lithologies[i].rhof
self.rhos[i] = self.lithologies[i].rhos
self.DeltaS[i] = self.lithologies[i].DeltaS
[docs] def bulkProperties(self, P=None, T=None):
"""
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
-------
dict
The bulk alpha, CP and rho for the mantle at the given P and T, labelled
as such.
"""
F = _np.zeros(self.number_lithologies)
if T is not None:
for i in range(self.number_lithologies):
F[i] = self.lithologies[i].F(P, T)
alpha = sum(self.proportions * self.alphaf * F + self.proportions * self.alphas * (1 - F))
CP = sum(self.proportions * self.CP)
rho = sum(self.proportions * self.rhof * F + self.proportions * self.rhos * (1 - F))
return {'alpha': alpha, 'CP': CP, 'rho': rho}
[docs] def solidusIntersection(self, Tp):
"""
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
-------
numpy.array
The pressure of solidus intersection of each lithology.
"""
intersect = _np.zeros(self.number_lithologies)
for i in range(self.number_lithologies):
def f_solve(P):
return self.lithologies[i].TSolidus(P) - self.adiabat(P, Tp)
# Check that the lithology actually melts:
if self.lithologies[i].TSolidus(3.0) is _np.inf:
intersect[i] = _np.nan
# Check there is actually some of the lithology:
elif self.proportions[i] > 0:
try:
result = _root_scalar(f_solve, x0=3.0, x1=4.0)
if result.converged is True:
intersect[i] = result.root
else:
intersect[i] = _np.nan
except Exception:
intersect[i] = _np.nan
else:
intersect[i] = _np.nan
return intersect
[docs] def solidusIntersectionIsobaric(self, P):
"""
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
-------
numpy.Array
The temperature of solidus intersection of each lithology.
"""
intersect = _np.zeros(self.number_lithologies)
for i in range(self.number_lithologies):
# Check whether there is actually any of the lithology present:
if self.proportions[i] > 0:
intersect[i] = self.lithologies[i].TSolidus(P)
else:
intersect[i] = _np.nan
return intersect
[docs] def adiabat(self, P, Tp):
"""
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.
"""
bulk_props = self.bulkProperties(P)
T = ((Tp + 273.15) * _np.exp(bulk_props['alpha']
/ (bulk_props['rho'] * bulk_props['CP']) * P) - 273.15)
return T
[docs] def F(self, P, T):
"""
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
-------
numpy.Array
Array containing the melt fraction of each lithology.
"""
F = _np.zeros(self.number_lithologies)
for i in range(self.number_lithologies):
F[i] = self.lithologies[i].F(P, T)
return F
[docs] def dFdP(self, P, T, prevent_freezing=False, F_prev=None):
"""
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
-------
numpy.Array
Array of dFdP values for each lithology.
"""
dFdP = _np.zeros(self.number_lithologies)
dTdP = _np.zeros(self.number_lithologies)
dTdF = _np.zeros(self.number_lithologies)
F = _np.zeros(self.number_lithologies)
for i in range(self.number_lithologies):
dTdP[i] = self.lithologies[i].dTdP(P, T)
dTdF[i] = self.lithologies[i].dTdF(P, T)
F[i] = self.lithologies[i].F(P, T)
lithologies_melting = _np.where(((F > 0) & (F < 1)) & (self.proportions > 0))[0]
if _np.shape(lithologies_melting)[0] > 0:
for i in range(_np.shape(lithologies_melting)[0]):
not_key = [True] * self.number_lithologies
key = lithologies_melting[i]
not_key[key] = False
bulk_props = self.bulkProperties(P, T)
# Equation (26) from PM2001 to find dFdP of first lithology
top = (bulk_props['CP'] / (T + 273.15) * dTdP[key]
- bulk_props['alpha'] / bulk_props['rho']
+ sum(self.proportions[not_key] * self.DeltaS[not_key]
* (dTdP[key] - dTdP[not_key]) / dTdF[not_key]))
bottom = (self.proportions[key] * self.DeltaS[key]
+ sum(self.proportions[not_key] * self.DeltaS[not_key]
* dTdF[key] / dTdF[not_key])
+ bulk_props['CP'] / (T + 273.15) * dTdF[key])
dFdP[key] = - top / bottom
if prevent_freezing and F[i] < F_prev[i]:
dFdP[key] = 0.0
return dFdP
[docs] def adiabaticGradient(self, P, T):
"""
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
-------
float
The adiabatic gradient in C/GPa
"""
bulk_props = self.bulkProperties(P, T)
dTdP = bulk_props['alpha'] * (T + 273.15) / bulk_props['rho'] / bulk_props['CP']
return dTdP
[docs] def dTdP(self, P, T, dFdP=None):
"""
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
-------
float
The thermal gradient in the melting region at the P and T of interest.
"""
if dFdP is None:
dFdP = self.dFdP(P, T)
melting_lithologies = _np.where(dFdP != 0)[0]
if _np.shape(melting_lithologies)[0] > 0:
key = _np.argmax(_np.abs(dFdP[dFdP != 0]))
key = melting_lithologies[key]
dTdP = self.lithologies[key].dTdP(P, T) + self.lithologies[key].dTdF(P, T) * dFdP[key]
else:
dTdP = self.adiabaticGradient(P, T)
return dTdP
[docs] def isobaricMelt(self, Tstart, P, dT=0.1):
"""
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
-------
float
Temperature of mantle following the melting step.
"""
solidus_intersection = self.solidusIntersectionIsobaric(P)
# Calculate the entropy lost associated with cooling solid material to the
# solidus temperature
solT = _np.nanmin(solidus_intersection)
bulk_props = self.bulkProperties(P, solT)
DeltaS_cool = - bulk_props['CP'] * _np.log((solT + 273.15) / (Tstart + 273.15))
DeltaS_melt = 0
T = solT + dT
while DeltaS_melt < DeltaS_cool and T < Tstart:
bulk_props = self.bulkProperties(P, T)
DeltaS_melt = (_np.sum(self.F(P, T) * self.proportions * self.DeltaS) +
- bulk_props['CP'] * _np.log((solT + 273.15) / (T + 273.15)))
T = T + dT
return T
[docs] def adiabaticMelt(self, Tp, Pstart=None, Pend=0.01, dP=-0.004, steps=None, ReportSSS=True,
adjust_pressure=True, prevent_freezing=True, warn_prevent_freezing=True):
"""
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
-------
pyMelt.meltingColumn
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.
"""
solidus_intersect = self.solidusIntersection(Tp)
if Pstart is None:
if all(_np.isnan(solidus_intersect)) is True:
raise InputError("No solidus intersection found. To model adiabatic "
"decompression of solid mantle set a starting pressure using "
"Pstart.")
Pstart = _np.nanmax(solidus_intersect) + 1e-5
adjust_pressure = False
if steps is None:
if dP >= 0:
raise InputError("dP should be less than zero.")
P = _np.arange(Pstart, Pend, dP)
steps = len(P)
else:
P = _np.linspace(Pstart, Pend, steps)
dP = (Pend - Pstart) / (steps - 1)
T = _np.zeros(steps)
T[0] = self.adiabat(Pstart, Tp)
if (adjust_pressure is True
and T[0] < _np.nanmin(self.solidusIntersectionIsobaric(Pstart))
and all(_np.isnan(solidus_intersect)) is False):
p_intersect = _np.nanmax(self.solidusIntersection(Tp))
diff = _np.abs(P - p_intersect)
adjustment = P[_np.argmin(diff)] - p_intersect
if adjustment < 0:
adjustment += dP
P -= adjustment
if P[-1] < 0:
P[-1] = 0
F = _np.zeros([steps, self.number_lithologies])
if T[0] > _np.nanmin(self.solidusIntersectionIsobaric(Pstart)):
if ReportSSS is True:
_warn("Super solidus start")
T[0] = self.isobaricMelt(T[0], Pstart)
for i in range(steps):
if i == 0:
F[i] = self.F(P[0], T[0])
else:
k1 = self.dFdP(P[i - 1], T[i - 1], prevent_freezing, F[i - 1])
j1 = self.dTdP(P[i - 1], T[i - 1], k1)
k2 = self.dFdP(P[i - 1] + dP / 2, T[i - 1] + dP / 2 * j1, prevent_freezing,
F[i - 1])
j2 = self.dTdP(P[i - 1] + dP / 2, T[i - 1] + dP / 2 * j1, k2)
k3 = self.dFdP(P[i - 1] + dP / 2, T[i - 1] + dP / 2 * j2, prevent_freezing,
F[i - 1])
j3 = self.dTdP(P[i - 1] + dP / 2, T[i - 1] + dP / 2 * j2, k3)
k4 = self.dFdP(P[i], T[i - 1] + dP * j3, prevent_freezing, F[i - 1])
j4 = self.dTdP(P[i], T[i - 1] + dP * j3, k4)
T[i] = T[i - 1] + dP / 6 * (j1 + 2 * j2 + 2 * j3 + j4)
F[i] = self.F(P[i], T[i])
if prevent_freezing is True:
for j in range(self.number_lithologies):
if F[i, j] < F[i - 1, j]:
F[i, j] = F[i - 1, j]
if warn_prevent_freezing is True:
_warn("Freezing prevented.")
results = _pd.DataFrame(F, columns=self.names)
results['P'] = P
results['T'] = T
return _meltingColumn(results, self, Tp)