"""
==================
McKenzie and Bickle (1988)
==================
Implementation of the anhydrous melting model presented by McKenzie and Bickle (1988).
"""
from pyMelt.lithology_classes import lithology as _lithology
import numpy as _np
from scipy.optimize import fsolve as _fsolve
from scipy.special import expit as _expit
from scipy.misc import derivative as _derivative
from scipy.optimize import root_scalar as _root_scalar
[docs]class lherzolite(_lithology):
"""
Implementation of the McKenzie and Bickle (1988) garnet peridotite melting model.
As this parameterisation provides pressure as a function of solidus temperature,
scipt.optimize.fsolve and scipy.special.expit are required to numerically find
solidus temperature as a function of pressure.
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 object, with
values:
- A1: Parameter used in solidus definition.
- A2: Parameter used in solidus definition.
- A3: Parameter used in solidus definition.
- A4: Parameter used in solidus definition.
- B1: Parameter used in liquidus definition.
- B2: Parameter used in liquidus definition.
- B3: Parameter used in liquidus definition.
- B4: Parameter used in liquidus definition.
- a0: Parameter used in calculating melt fraction.
- a1: Parameter used in calculating melt fraction.
The thermal expansivities, the heat capacity, the densities, and the entropy of fusion may
also be changed during object initialisation.
Parameters
----------
CP : float, default: 1120.0
The heat capacity (J K-1 kg-1)
alphas : float, default: 40
The thermal expansivity of the solid (1e-6 K-1)
alphaf : float, default: 68
The thermal expansivity of the melt (1e-6 K-1)
rhos : float, default: 3.3
The density of the solid (kg m-3)
rhof : float, default: 2.8
The density of the melt (kg m-3)
DeltaS : float, default: 250
The entropy of fusion J K-1 kg-1
parameters : dict, default: parameters from McKenzie and Bickle (1988)
The model parameters described above
"""
def __init__(self,
CP=1120.0,
alphas=40,
alphaf=68,
rhos=3.3,
rhof=2.8,
DeltaS=250,
parameters={'A1': 1100.0,
'A2': 136.0,
'A3': 4.968e-4,
'A4': 1.2e-2,
'B1': 1736.0,
'B2': 4.343,
'B3': 180.0,
'B4': 2.2169,
'a0': 0.4256,
'a1': 2.988
}
):
self.DeltaS = DeltaS
self.CP = CP
self.alphas = alphas
self.alphaf = alphaf
self.rhos = rhos
self.rhof = rhof
self.parameters = parameters
[docs] def TSolidus(self, P):
"""
Calculates the solidus temperature, at a given pressure, using Equation 18 of
McKenzie and Bickle (1988). Requires scipy.optimize.fsolve to solve for pressure.
As fsolve can only solve for a single pressure TSolidus first detects whether P is
given as a list, tuple, numpy.ndarray (e.g. numpy.linspace), or scalar prior to
calculating the solidus pressure.
Parameters
----------
P : float
Pressure (GPa)
Returns
-------
float or numpy.Array
Solidus temperature (degC).
"""
TSolidus_initial_guess = 1300.0
def func(TSolidus):
_func = ((TSolidus - self.parameters['A1']) / self.parameters['A2']
+ self.parameters['A3'] * _np.exp(self.parameters['A4']
* (TSolidus - self.parameters['A1']))) - Pressure
return _func
if isinstance(P, (list, tuple, _np.ndarray)) is True:
TSolidus_solution = []
for i in range(len(P)):
Pressure = P[i]
TSolidus_solution.append(_fsolve(func, TSolidus_initial_guess, maxfev=50000)[0])
TSolidus_solution = _np.array(TSolidus_solution)
else:
Pressure = P
TSolidus_solution = _fsolve(func, TSolidus_initial_guess, maxfev=50000)[0]
return TSolidus_solution
[docs] def TLiquidus(self, P):
"""
Calculates the liquidus temperature, at a given pressure, using Equation 19 of
McKenzie and Bickle (1988).
Parameters
----------
P : float
Pressure (GPa)
Returns
-------
float or numpy.Array
Solidus temperature (degC).
"""
TLiquidus = (self.parameters['B1'] + self.parameters['B2'] * P + self.parameters['B3']
* _np.arctan(P / self.parameters['B4']))
return TLiquidus
def _dTdPSolidus(self, P):
"""
Returns the solidus temperature gradient at any given pressure. Requires
scipy.optimize.fsolve to solve for pressure, and scipy.special.expit to
avoid exponential overflows that can occur with numpy.exp.
Parameters
----------
P : float
Pressure (GPa).
Returns
-------
float
Solidus temperaure gradient (degC/GPa)
"""
TSolidus = self.TSolidus(P)
dTdPSolidus = (1 / (1 / self.parameters['A2'] + self.parameters['A3']
* self.parameters['A4'] * _expit(self.parameters['A4']
* (TSolidus - self.parameters['A1']))))
return dTdPSolidus
def _dTdPLiquidus(self, P):
"""
Returns the liquidus temperature gradient at any given pressure.
Parameters
----------
P : float
Pressure (GPa).
Returns
-------
float
Solidus temperaure gradient (degC/GPa)
"""
dTdPLiquidus = (self.parameters['B2'] + self.parameters['B3']
/ (1 + (P / self.parameters['B4'])**2))
return dTdPLiquidus
def _RescaledT(self, T, P):
"""
Calculates the rescaled temperature defined by Equation 20 of
McKenzie and Bickle (1988).
Parameters
----------
T : float
Temperature (degC).
P : float
Pressure (GPa).
Returns
-------
float
Rescaled Temperature (dimensionless).
"""
TSolidus = self.TSolidus(P)
TLiquidus = self.TLiquidus(P)
RescaledT = (T - (TSolidus + TLiquidus) / 2) / (TLiquidus - TSolidus)
return RescaledT
[docs] def F(self, P, T):
"""
Wrapper for the melt fraction functions, defined by Equations 21 and 22
of McKenzie and Bickle (1988). If T and P are below the solidus, returns 0,
if they are above the liquidus, returns 1. Otherwise determines F.
Parameters
----------
P : float
Pressure (GPa).
T : float
Temperature (degC).
Returns
-------
float
Melt fraction.
"""
TSolidus = self.TSolidus(P)
TLiquidus = self.TLiquidus(P)
RescaledT = self._RescaledT(T, P)
if T > TLiquidus:
F = 1.0
elif T < TSolidus:
F = 0.0
else:
F = (RescaledT + (RescaledT**2 - 0.25) * (self.parameters['a0']
+ self.parameters['a1'] * RescaledT) + 0.5)
return F
[docs] def dTdF(self, P, T, **kwargs):
"""
Calculates dT/dF(const. P) numerically.
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)
F = self.F(P, T, **kwargs)
if F == 0:
dTdF = _np.inf
elif F < 1:
dTdF = 1.0 / (_derivative(_to_diff, T, dx=0.001, args=(P, kwargs)))
else:
dTdF = _np.inf
return dTdF
[docs] def dTdP(self, P, T, **kwargs):
"""
Calculates dT/dP(const. F) numerically.
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