"""
=================
Lithology classes
=================
The lithology classes module provides the `lithology` class for implementing melting models, and
the `hydrous_lithology` class for converting an anhydrous lithology to a hydrous lithology.
"""
from scipy.misc import derivative as _derivative
from scipy.optimize import root_scalar as _root_scalar
import numpy as _np
from pyMelt.core import ConvergenceError
from warnings import warn as _warn
# Default constant values taken from Katz et al., 2003:
default_properties = {'CP': 1000.0, # Heat capacity in J Kg-1 K-1
'alphas': 40.0, # Thermal expansivity of the solid (K-1).
'alphaf': 68.0, # Thermal expansivity of the melt (K-1).
'rhos': 3.3, # Density of the solid (g cm-3).
'rhof': 2.9, # Density of the melt (g cm-3).
'DeltaS': 300.0, # Entropy of fusion. (J kg-1 K-1).
}
[docs]class lithology(object):
"""
Lithology base class. This class contains all the parameters and methods required to calculate
the melting behaviour of a single mantle component.
Parameters
----------
CP: float
Heat capacity (J Kg-1 K-1). Default is 1000.0.
alphas: float
Thermal expansivity of the solid (K-1). Default is 40.0.
alphaf: float
Thermal expansivity of the melt (K-1). Default is 68.0.
rhos: float
Density of the solid (g cm-3). Default is 3.3.
rhof: float
Density of the melt (g cm-3). Default is 2.9.
DeltaS: float
Entropy of fusion. (J kg-1 K-1). Default is 300.0.
parameters: dict
A dictionary of the parameters required by the lithology.
"""
def __init__(self,
CP=default_properties['CP'],
alphas=default_properties['alphas'],
alphaf=default_properties['alphaf'],
rhos=default_properties['rhos'],
rhof=default_properties['rhof'],
DeltaS=default_properties['DeltaS'],
parameters={}):
self.CP = CP
self.alphas = alphas
self.alphaf = alphaf
self.rhos = rhos
self.rhof = rhof
self.DeltaS = DeltaS
self.parameters = parameters
[docs] def TSolidus(self, P):
"""
Returns the temperature of the solidus at any given pressure.
Parameters
----------
P: float
Pressure (GPa).
Returns
-------
None
The default lithology never melts.
"""
return None
[docs] def TLiquidus(self, P):
"""
Returns the temperature of the liquidus at any given pressure.
Parameters
----------
P: float
Pressure (GPa).
Returns
-------
None
The default lithology never melts.
"""
return None
[docs] def F(self, P, T):
"""
Returns the melt fraction at any given pressure and temperature.
Parameters
----------
P: float
Pressure (GPa).
T: float
Temperature (degC).
Returns
-------
float
Melt fraction.
"""
return 0.0
[docs] def dTdF(self, P, T, **kwargs):
"""
Calculates dT/dF(const. P) numerically. For faster and more accurate results redefine
this method.
Parameters
----------
P: float
Pressure (GPa)
T: float
Temperature (degC)
Returns
-------
float
dT/dF(const. P) (K).
"""
def _to_diff(T, P, kwargs={}):
return self.F(P, T, **kwargs)
return 1.0 / (_derivative(_to_diff, T, dx=0.001, args=(P, kwargs)))
[docs] def dTdP(self, P, T, **kwargs):
"""
Calculates dT/dP(const. F) numerically. For faster and more accurate results redefine
this method.
Parameters
----------
P: float
Pressure (GPa).
T: float
Temperature (degC).
Returns
-------
float
dT/dP(const. F) (K GPa-1).
"""
# This finds the temperature at a given pressure for which the melt fraction is eqal to the
# value specified. Used for calculating dT/dP (at const. F).
def _to_diff_dTdP(P, F, T, kwargs={}):
t = _root_scalar(_hold_constant_F, x0=T, x1=T + 10, args=(P, F, kwargs)).root
return t
# This method is used to find the P-T curve at which F remains constant, for the
# calculation of dT/dP (at const. F).
def _hold_constant_F(t, P, F, kwargs={}):
return self.F(P, t, **kwargs) - F
if self.F(P, T, **kwargs) == 0:
dTdP = self.alphas / self.rhos / self.CP
elif self.F(P, T, **kwargs) == 1:
dTdP = self.alphas / self.rhos / self.CP
else:
F = self.F(P, T, **kwargs)
dTdP = _derivative(_to_diff_dTdP, P, dx=0.001, args=(F, T, kwargs))
return dTdP
[docs]class hydrousLithology(object):
r"""
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):
.. math::
\Delta T(X_{H2O}) = K X^\gamma_{H2O}
Where :math:`X_{H2O}` is the water content in the melt.
The water content at saturation is controlled by eqn. 17 of Katz et al. (2003):
.. math::
X^{sat}_{H2O} = \chi_1 P^\lambda + \chi_2 P
"""
def __init__(self, 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):
self.lithology = lithology
for m in dir(lithology):
if '__' not in m and m not in ['TSolidus', 'TLiquidus', 'F', 'dTdF', 'dTdP']:
f = lithology.__getattribute__(m)
if callable(f):
self.__setattr__(m, f.__func__.__get__(self, self.__class__))
else:
self.__setattr__(m, f)
self.F_function = self.lithology.F.__func__.__get__(self, self.__class__)
self.H2O = H2O
self.K = K
self.gamma = gamma
self.D = D
self.chi1 = chi1
self.chi2 = chi2
self.l = l
self.continuous = continuous
self.phi = phi
self.threshold_F = threshold_F
[docs] def H2O_saturation(self, P):
"""
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
-------
float
The H2O content of the melt at H2O-saturation, in wt%.
"""
return self.chi1 * P ** self.l + self.chi2 * P
[docs] def melt_H2O(self, F, P=None):
"""
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
-------
float
The H2O content of the melt, in wt%.
"""
if self.continuous is False:
Cl = self.H2O / ((1.0 - F) * self.D + F)
else:
Cl = (self.H2O / ((1 - self.phi) * self.D + self.phi)
* (1 - F)**((1 - self.phi) * (1 - self.D) / ((1 - self.phi) * self.D)))
# Test for H2O saturation:
if P is not None:
H2Osat = self.H2O_saturation(P)
if H2Osat < Cl:
if self.continuous:
_warn("Fluid saturated during continuous melting- you may get bizarre "
"behaviour")
Cl = H2Osat
return Cl
[docs] def TSolidus(self, P, F=0.0):
"""
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
-------
float
Solidus temperature (degC).
"""
TSolidus = self.lithology.TSolidus(P) - self.K * self.melt_H2O(F, P)**self.gamma
return TSolidus
[docs] def TLiquidus(self, P, **kwargs):
"""
Returns the temperature of the liquidus at any given pressure.
Parameters
----------
P: float
Pressure (GPa).
Returns
-------
float
Liquidus temperature (degC).
"""
return self.lithology.TLiquidus(P)
[docs] def F(self, P, T, **kwargs):
"""
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
-------
float
The melt fraction.
"""
if self.threshold_F is not None:
anhydrous_F = self.lithology.F(P, T, **kwargs)
if anhydrous_F > self.threshold_F:
return anhydrous_F
F = _root_scalar(self._f_to_solve, bracket=[0, 1], args=(P, T))
if F.flag != 'converged':
raise ConvergenceError('The melt fraction calculation did not converge.')
return _np.nan
else:
return F.root
[docs] def dTdF(self, P, T, **kwargs):
"""
Calculates dT/dF(const. P) for the hydrated lithology numerically.
Parameters
----------
P: float
Pressure (GPa)
T: float
Temperature (degC)
Returns
-------
float
dT/dF(const. P) (K).
"""
if self.F(P, T, **kwargs) == 0 or self.F(P, T, **kwargs) == 1:
return _np.inf
return 1.0 / (_derivative(self._to_diff_dTdF, T, dx=0.001, args=(P, kwargs)))
[docs] def dTdP(self, P, T, **kwargs):
"""
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
-------
float
dT/dP(const. F) (K GPa-1).
"""
if self.F(P, T, **kwargs) == 0:
dTdP = self.alphas / self.rhos / self.CP
elif self.F(P, T, **kwargs) == 1:
dTdP = self.alphas / self.rhos / self.CP
else:
F = self.F(P, T, **kwargs)
dTdP = _derivative(self._to_diff_dTdP, P, dx=0.001, args=(F, T, kwargs))
return dTdP
# The following methods support the numerical differentiation and root-finding methods used
# above. The reason they are defined as methods, and do not sit within their "parent" methods
# is to assist with module testing and troubleshooting.
# This method rearranges the arguments so that F can be numerically differentiated.
def _to_diff_dTdF(self, T, P, kwargs={}):
return self.F(P, T, **kwargs)
# This method is used to find the P-T curve at which F remains constant, for the calculation of
# dT/dP (at const. F).
def _hold_constant_F(self, t, P, F, kwargs={}):
return self.F(P, t, **kwargs) - F
# This finds the temperature at a given pressure for which the melt fraction is eqal to the
# value specified. Used for calculating dT/dP (at const. F).
def _to_diff_dTdP(self, P, F, T, kwargs={}):
t = _root_scalar(self._hold_constant_F, x0=T, x1=T + 10, args=(P, F, kwargs)).root
return t
# This provides the misfit function when solving for F (Eqn. 19 of Katz et al., 2003)
def _f_to_solve(self, x, P, T):
fcalc = self.F_function(P, T, F=x)
misfit = fcalc - x
return misfit