Source code for pyMelt.lithologies.shorttle

"""
======================
Shorttle et al. (2014)
======================

Implementation of the new lithologies in Shorttle et al. (2014).

"""

from pyMelt.lithology_classes import lithology as _lithology

import numpy as _np


[docs]class kg1(_lithology): """ Implementation of the KG1 parameterisation by Shorttle et al. (2014). 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: Constant used in solidus expression. - A2: Constant used in solidus expression. - A3: Constant used in solidus expression. - B1: Constant used in cpx-out expression. - B2: Constant used in cpx-out expression. - B3: Constant used in cpx-out expression. - C1: Constant used in liquidus expression. - C2: Constant used in liquidus expression. - C3: Constant used in liquidus expression. - a: Constant used in cpx-present melt fraction expression. - b: Constant used in cpx-present melt fraction expression. - c: Constant used in cpx-absent melt fraction expression. - d: Constant used in cpx-absent melt fraction expression. - alpha: Exponent used in the cpx-present melt fraction expression. - beta: Exponent used in the cpx-absent melt fraction expression. The thermal expansivities, the heat capacity, the densities, and the entropy of fusion may also be changed during object initialisation. Parameters ---------- CP : float, default: 1140.0 The heat capacity (J K-1 kg-1) alphas : float, default: 30.0 The thermal expansivity of the solid (1e-6 K-1) alphaf : float, default: 68.0 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.9 The density of the melt (kg m-3) DeltaS : float, default: 380.0 The entropy of fusion J K-1 kg-1 parameters : dict, default: parameters from Shorttle et al. (2014) The model parameters described above """ def __init__(self, CP=1140.0, alphas=30.0, alphaf=68.0, rhos=3.3, rhof=2.9, DeltaS=380.0, parameters={'A1': 1095.4, 'A2': 124.1, 'A3': - 4.7, 'B1': 1179.6, 'B2': 157.2, 'B3': - 11.1, 'C1': 1780.0, 'C2': 45.0, 'C3': - 2.0, 'a': 0.3187864, 'b': 0.4154, 'c': 0.7341864, 'd': 0.2658136, 'alpha': 2.0, 'beta': 1.5 } ): 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, **kwargs): """ Returns solidus temperature at a given pressure. T = A1 + A2*P + A3*P**2. Parameters ---------- P : float, or list like Pressure in GPa. Returns ------- float, or list like Solidus Temperature in degC. """ T = self.parameters['A1'] + self.parameters['A2'] * P + self.parameters['A3'] * P**2 return T
[docs] def TLiquidus(self, P, **kwargs): """ Returns liquidus temperature at a given pressure. T = C1 + C2*P + C3*P**2. Parameters ---------- P : float, or list like Pressure in GPa. Returns ------- float, or list like Liquidus temperature in degC. """ T = self.parameters['C1'] + self.parameters['C2'] * P + self.parameters['C3'] * P**2 return T
[docs] def dTdP(self, P, T, **kwargs): """ Returns dT/dP (constant F) at a given pressure and temperature. Parameters ---------- P : float Pressure in GPa. T : float Temperature in degC. Returns ------- float dT/dP (constant F) in K GPa-1. """ Tsol = self.TSolidus(P, **kwargs) Tliq = self.TLiquidus(P, **kwargs) Tcpx = self._TCpxOut(P, **kwargs) if T < Tsol: dTdP = self.alphas / self.rhos / self.CP if T < Tcpx: dTdP = (-(-(T - Tsol) / (Tcpx - Tsol) * (self.parameters['B2'] + 2 * self.parameters['B3'] * P - self.parameters['A2'] - self.parameters['A3'] * 2 * P) - self.parameters['A2'] - self.parameters['A3'] * 2 * P)) elif T < Tliq: dTdP = (-(-(T - Tcpx) / (Tliq - Tcpx) * (self.parameters['C2'] + self.parameters['C3'] * 2 * P - self.parameters['B2'] - self.parameters['B3'] * 2 * P) - self.parameters['B2'] - 2 * self.parameters['B3'] * P)) else: dTdP = self.alphas / self.rhos / self.CP return dTdP
[docs] def dTdF(self, P, T, **kwargs): """ Returns dT/dF (constant P) at a given pressure and temperature. If below the solidus, or above the liquidus, _np.inf is returned. Parameters ---------- P : float Pressure in GPa. T : float Temperature in degC. Returns ------- float dT/dF (constant P) in K. """ Tsol = self.TSolidus(P, **kwargs) Tliq = self.TLiquidus(P, **kwargs) Tcpx = self._TCpxOut(P, **kwargs) if T < Tsol: dTdF = _np.inf elif T < Tcpx: dTdF = ((Tcpx - Tsol) / (self.parameters['b'] + self.parameters['a'] * self.parameters['alpha'] * ((T - Tsol) / (Tcpx - Tsol))**(self.parameters['alpha'] - 1))) elif T < Tliq: dTdF = ((Tliq - Tcpx) / (self.parameters['d'] * self.parameters['beta'] * ((T - Tcpx) / (Tliq - Tcpx))**(self.parameters['beta'] - 1))) else: dTdF = _np.inf return dTdF
[docs] def F(self, P, T, **kwargs): """ Returns melt fraction at a given pressure and temperature. If below the solidus, returns 0. If above the liquidus, returns 1. Prior to cpx exhaustion: F = a*(Tr)**alpha _ b*Tr Tr = (T-Tsol)/(Tliq-Tsol) After cpx exhaustion: F = d*(Tr)**beta + c Tr = (T-Tcpx)/(Tliq-Tcpx) Parameters ---------- P : float Pressure in GPa. T : float Temperature in degC. Returns ------- float Melt fraction between 0 and 1. """ Tsol = self.TSolidus(P, **kwargs) Tliq = self.TLiquidus(P, **kwargs) Tcpx = self._TCpxOut(P, **kwargs) if T < Tsol: F = 0.0 elif T > Tliq: F = 1.0 elif T < Tcpx: Tf = (T - Tsol) / (Tcpx - Tsol) F = self.parameters['a'] * Tf**self.parameters['alpha'] + self.parameters['b'] * Tf else: Tf = (T - Tcpx) / (Tliq - Tcpx) F = self.parameters['d'] * Tf**self.parameters['beta'] + self.parameters['c'] return F
def _TCpxOut(self, P, **kwargs): """ Returns the temperature of cpx-exhaustion at a given pressure. T = B1 + B2*P + B3*P**2. Parameters ---------- P :float, or list like Pressure in GPa. Returns ------- float, or list like Cpx-exhaustion temperature in degC. """ T = self.parameters['B1'] + self.parameters['B2'] * P + self.parameters['B3'] * P**2 return T
[docs]class harzburgite(_lithology): """ Material that does not melt, i.e. Harzburgite in Shorttle et al. (2014) and Matthews et al. (2016). Default constants as defined by Shorttle et al. (2014). The thermal expansivities, the heat capacity, the densities, and the entropy of fusion may be changed during object initialisation. Parameters ---------- CP : float, default: 1000.0 The heat capacity (J K-1 kg-1) alphas : float, default: 30.0 The thermal expansivity of the solid (1e-6 K-1) alphaf : float, default: 30.0 Melt thermal expansivity, not used, here for consistency. rhos : float, default: 3.25 The density of the solid (kg m-3) rhof : float, default: 2.9 Melt density, not used, here for consistency. DeltaS : float, default: 300 The entropy of fusion, not used, here for consistency. parameters : dict, default: {} This model does not use any parameters, here for consistency. """ def __init__(self, CP=1000.0, alphas=30.0, alphaf=30.0, rhos=3.25, rhof=2.9, DeltaS=300, parameters={} ): self.CP = CP self.alphas = alphas self.alphaf = alphaf self.rhos = rhos self.rhof = rhof self.DeltaS = DeltaS self.parameters = parameters
[docs] def F(self, P, T, **kwargs): """ Melt Fraction. Returns 0.0. Parameters ---------- P : float Pressure. There to maintain consistancy within lithology definitions. T : float Temperature. There to maintain consistancy within lithology definitions. Returns ------- float Melt fraction will always be 0.0. """ return 0.0
[docs] def dTdF(self, P, T, **kwargs): """ dTdF(constP). Returns _np.inf. Parameters ---------- P : float Pressure. There to maintain consistancy within lithology definitions. T : float Temperature. There to maintain consistancy within lithology definitions. Returns ------- numpy.inf The value will always be infinite. """ return _np.inf
[docs] def dTdP(self, P, T, **kwargs): """ dTdP(constF). Returns 0.0. Parameters ---------- P : float Pressure. There to maintain consistancy within lithology definitions. T : float Temperature. There to maintain consistancy within lithology definitions. Returns ------- float The value will always be 0.0. """ dTdP = self.alphas / self.rhos / self.CP return dTdP
[docs] def TSolidus(self, P, **kwargs): """ Solidus temperature. Returns _np.inf. Parameters ---------- P : float Pressure. There to maintain consistancy within lithology definitions. Returns ------- numpy.inf The value will always be infinite. """ if isinstance(P, list) or isinstance(P, _np.ndarray): return _np.array([_np.inf] * len(P)) else: return _np.inf
[docs] def TLiquidus(self, P, **kwargs): """ Liquidus temperature. Returns _np.inf Parameters ---------- P : float Pressure. There to maintain consistancy within lithology definitions. Returns ------- numpy.inf The value will always be infinite. """ if isinstance(P, list) or isinstance(P, _np.ndarray): return _np.array([_np.inf] * len(P)) else: return _np.inf