Source code for pyMelt.lithologies.katz

"""
==================
Katz et al. (2003)
==================

Implementation of the anhydrous melting model presented by Katz et al. (2003).

"""

from pyMelt.lithology_classes import lithology as _lithology
from pyMelt.lithology_classes import default_properties as _default_properties

import numpy as _np


[docs]class lherzolite(_lithology): """ Implementation of the Katz et al. (2003) anhydrous lherzolite melting model. 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. - Mcpx: Mass fraction of cpx in the source. Controls the transition to low-productivity harzburgite-type melting. - A1: Parameter used to define solidus. - A2: Parameter used to define solidus. - A3: Parameter used to define solidus. - B1: Parameter used to define lherzolite-liquidus. - B2: Parameter used to define lherzolite-liquidus. - B3: Parameter used to define lherzolite-liquidus. - C1: Parameter used to define liquidus. - C2: Parameter used to define liquidus. - C3: Parameter used to define liquidus. - beta1: Parameter used to calculate melt fraction during cpx-present melting. - beta2: Parameter used to calculate melt fraction during cpx-absent melting. - r1: Parameter used to define cpx reaction coefficient. - r2: Parameter used to define cpx reaction coefficient. The thermal expansivities, the heat capacity, the densities, and the entropy of fusion may also be changed during object initialisation. Parameters ---------- CP : float, default: pyMelt.lithology_class.default_properties['CP'] The heat capacity (J K-1 kg-1) alphas : float, default: pyMelt.lithology_class.default_properties['alphas'] The thermal expansivity of the solid (1e-6 K-1) alphaf : float, default: pyMelt.lithology_class.default_properties['alphaf'] The thermal expansivity of the melt (1e-6 K-1) rhos : float, default: pyMelt.lithology_class.default_properties['rhos'] The density of the solid (kg m-3) rhof : float, default: pyMelt.lithology_class.default_properties['rhof'] The density of the melt (kg m-3) DeltaS : float, default: pyMelt.lithology_class.default_properties['DeltaS'] The entropy of fusion J K-1 kg-1 parameters : dict, default: parameters from Katz et al. (2003) The model parameters described above """ 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={'Mcpx': 0.17, 'A1': 1085.70, 'A2': 132.9, 'A3': - 5.1, 'B1': 1475.0, 'B2': 80.0, 'B3': - 3.2, 'C1': 1780.0, 'C2': 45.0, 'C3': - 2.0, 'beta1': 1.5, 'beta2': 1.5, 'r1': 0.5, 'r2': 0.08 } ): 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 the temperature of the solidus at any given pressure. Eqn(4). Parameters ---------- P : float Pressure (GPa). Returns ------- float Solidus temperature (degC). """ TSolidus = (self.parameters['A1'] + self.parameters['A2'] * P + self.parameters['A3'] * (P**2)) return TSolidus
[docs] def TLiquidus(self, P, **kwargs): """ Returns the temperature of the liquidus at any given pressure. Eqn(10). Parameters ---------- P : float Pressure (GPa). Returns ------- float Liquidus temperature (degC). """ TLiquidus = (self.parameters['C1'] + self.parameters['C2'] * P + self.parameters['C3'] * (P**2)) return TLiquidus
[docs] def F(self, P, T, **kwargs): """ Wrapper for the melt fraction functions. If T and P are below the solidus, returns 0, if they are above the liquidus, returns 1. If below the temperature of cpx-exhaustion, calls the Fcpx function, otherwise calls the Fopx function. Parameters ---------- P : float Pressure (GPa). T : float Temperature (degC). Returns ------- float Melt fraction. """ if T > self.TLiquidus(P, **kwargs): F = 1.0 elif T < self.TSolidus(P, **kwargs): F = 0.0 elif T < self._TcpxOut(P, **kwargs): F = self._Fcpx(T, P, **kwargs) else: F = self._Fopx(T, P, **kwargs) return F
[docs] def dTdF(self, P, T, **kwargs): """ Calculates dT/dF(const. P). First calculates the melt fraction. If F is zero, returns _np.inf. If F is 1, returns _np.inf. Otherwise uses the appropriate expressions for cpx present or absent melting. Parameters ---------- P : float Pressure (GPa) T : float Temperature (degC) Returns ------- float dT/dF(const. P) (K). """ F = self.F(P, T, **kwargs) if F == 0: # If no melt fraction the derivative is zero. Prevents division by zero. dTdF = _np.inf elif F < self._FcpxOut(P, **kwargs): dTdF = (((1 / self.parameters['beta1'])) * (self._TLherzLiquidus(P, **kwargs) - self.TSolidus(P, **kwargs)) * (F**((1 - self.parameters['beta1']) / self.parameters['beta1']))) elif F < 1.0: dTdF = (((1 / self.parameters['beta2'])) * (self.TLiquidus(P, **kwargs) - self._TcpxOut(P, **kwargs)) * (F**((1 - self.parameters['beta2']) / self.parameters['beta2']))) else: dTdF = _np.inf return dTdF
[docs] def dTdP(self, P, T, **kwargs): """ Calculates dT/dP(const. F). First calculates F, then chooses the appropriate expression for cpx present or absent melting. Parameters ---------- P : float Pressure (GPa). T : float Temperature (degC). Returns ------- float dT/dP(const. F) (K GPa-1). """ F = self.F(P, T, **kwargs) FcpxOut = self._FcpxOut(P, **kwargs) dTdPSolidus = self._dTdPSolidus(P, **kwargs) TLiquidus = self.TLiquidus(P, **kwargs) dTdPLiquidus = self._dTdPLiquidus(P, **kwargs) dTdPLherzLiquidus = self._dTdPLherzLiquidus(P, **kwargs) TcpxOut = self._dTdPcpxOut(P, **kwargs) dTdPcpxOut = self._dTdPcpxOut(P, **kwargs) FcpxOut = self._FcpxOut(P, **kwargs) dFdPcpxOut = self._dFdPcpxOut(P, **kwargs) if F == 0: dTdP = self.alphas / self.rhos / self.CP elif F < self._FcpxOut(P, **kwargs): dTdP = (((F**(1 / self.parameters['beta1'])) * (dTdPLherzLiquidus - dTdPSolidus)) + dTdPSolidus) elif F < 1.0: Trel = (T - TcpxOut) / (TLiquidus - TcpxOut) dTdP = ((TLiquidus - TcpxOut) / (1 - FcpxOut) * (1 / self.parameters['beta2']) * Trel**(1 - self.parameters['beta2']) * dFdPcpxOut * (Trel**self.parameters['beta2'] - 1) + dTdPcpxOut + Trel * (dTdPLiquidus - dTdPcpxOut)) else: dTdP = self.alphaf / self.rhof / self.CP return dTdP
def _TLherzLiquidus(self, P, **kwargs): """ Returns the temperature of the lherzolite liquidus at any given pressure. Eqn(5). This is the temperature at which the rock would be completely molten if cpx was remained present for the entirety of melting. Parameters ---------- P : float Pressure (GPa). Returns ------- Tlzliq: float Lherzolite liquidus temperature (degC). """ TLherzLiquidus = (self.parameters['B1'] + self.parameters['B2'] * P + self.parameters['B3'] * (P**2)) return TLherzLiquidus def _RescaledTcpx(self, T, P, **kwargs): """ Calculates the rescaled temperature during cpx-present melting. (Eq 3). Parameters ---------- T : float Temperature (degC). P : float Pressure (GPa). Returns ------- float Rescaled Temperature (dimensionless). """ TSolidus = self.TSolidus(P, **kwargs) TLherzLiquidus = self._TLherzLiquidus(P, **kwargs) RescaledTemperaturecpx = ((T - TSolidus) / (TLherzLiquidus - TSolidus)) return RescaledTemperaturecpx def _Fcpx(self, T, P, **kwargs): """ Melt fraction during cpx-present melting at the given P and T. Eq(2). Parameters ---------- T : float Temperature (degC). P : float Pressure (degC). Returns ------- float Melt fraction during cpx-present melting. """ RescaledT = self._RescaledTcpx(T, P, **kwargs) Fcpx = RescaledT**self.parameters['beta1'] return Fcpx def _RxnCoef(self, P, **kwargs): """ Reaction coefficient for cpx during melting at the specified pressure. Eq(7). Parameters ---------- P : float Pressure (GPa). Returns ------- float Reaction Coefficient. """ RxnCoef = self.parameters['r1'] + self.parameters['r2'] * P return RxnCoef def _FcpxOut(self, P, **kwargs): """ Calculates the melt fraction required to exhaust cpx from the residue at the given pressure. Eq(6). Parameters ---------- P : float Pressure (GPa) Returns ------- float Melt fraction at which cpx is exhausted. """ FcpxOut = self.parameters['Mcpx'] / self._RxnCoef(P, **kwargs) return FcpxOut def _TcpxOut(self, P, **kwargs): """ Calculates the temperature at which cpx will be exhausted during melting at the given pressure. Eq(9). Parameters ---------- P : float Pressure (GPa). Returns ------- float Temperature of cpx-exhaustion. """ TSolidus = self.TSolidus(P, **kwargs) TcpxOut = (((self._FcpxOut(P, **kwargs)**(1 / self.parameters['beta1']))) * (self._TLherzLiquidus(P, **kwargs) - TSolidus) + TSolidus) return TcpxOut def _RescaledTopx(self, T, P, **kwargs): """ Calculates the rescaled temperature during cpx-absent melting. Parameters ---------- T : float Temperature (degC). P : float Pressure (GPa). Returns ------- float Rescaled Temperature (dimensionless). """ TcpxOut = self._TcpxOut(P, **kwargs) RescaledTopx = ((T - TcpxOut) / (self.TLiquidus(P, **kwargs) - TcpxOut)) return RescaledTopx def _Fopx(self, T, P, **kwargs): """ Melt fraction during cpx-absent melting at the given P and T. Eq(8). Parameters ---------- T : float Temperature (degC). P : float Pressure (degC). Returns ------- float Melt fraction during cpx-absent melting. """ FcpxOut = self._FcpxOut(P, **kwargs) FopxDry = (FcpxOut + (1 - FcpxOut) * self._RescaledTopx(T, P, **kwargs)**self.parameters['beta2']) return FopxDry def _dTdPSolidus(self, P, **kwargs): return self.parameters['A2'] + 2 * self.parameters['A3'] * P def _dTdPLiquidus(self, P, **kwargs): return self.parameters['C2'] + 2 * self.parameters['C3'] * P def _dTdPLherzLiquidus(self, P, **kwargs): return self.parameters['B2'] + 2 * self.parameters['B3'] * P def _dTdPcpxOut(self, P, **kwargs): # term1 = ((self._TLherzLiquidus(P, **kwargs) - self.TSolidus(P, **kwargs)) # / self.parameters['beta1'] # * self._FcpxOut(P, **kwargs)**(1 / self.parameters['beta1'] - 1)) # term2 = (self._FcpxOut(P, **kwargs)**(1 / self.parameters['beta1']) # * (self._dTdPLherzLiquidus(P, **kwargs) - self._dTdPSolidus(P, **kwargs)) # + self._dTdPSolidus(P, **kwargs)) # return term1 * self._dFdPcpxOut(P, **kwargs) + term2 TSolidus = self.TSolidus(P, **kwargs) TLherzLiquidus = self._TLherzLiquidus(P, **kwargs) FcpxOut = self._FcpxOut(P, **kwargs) dTdPSolidus = self._dTdPSolidus(P, **kwargs) dTdPLherzLiquidus = self._dTdPLherzLiquidus(P, **kwargs) dFdPcpxOut = self._dFdPcpxOut(P, **kwargs) A = ((dFdPcpxOut * (1 / self.parameters['beta1']) * ((FcpxOut**((1 / self.parameters['beta1']) - 1)))) * (TLherzLiquidus - TSolidus)) B = ((FcpxOut**(1 / self.parameters['beta1'])) * (dTdPLherzLiquidus - dTdPSolidus)) + dTdPSolidus dTdPcpxOut = A + B return dTdPcpxOut def _dFdPcpxOut(self, P, **kwargs): # return - self.parameters['Mcpx'] / self._RxnCoef(P, **kwargs)**2 * self.parameters['r2'] RxnCoef = self._RxnCoef(P, **kwargs) FcpxOut = self._FcpxOut(P, **kwargs) dFdPcpxOut = ((-1) * FcpxOut * self.parameters['r2']) / RxnCoef return dFdPcpxOut