Source code for pyMelt.chemistry

"""
=========
Chemistry
=========

The chemistry module provides the base classes for defining chemical elements/species for
inclusion in pyMelt calculations, alongside default implementations.
"""

import numpy as _np
import pandas as _pd
from warnings import warn as _warn

default_methods = {
    "Rb": "continuous_instantaneous",
    "Ba": "continuous_instantaneous",
    "Th": "invmel",
    "U": "invmel",
    "Nb": "invmel",
    "Ta": "invmel",
    "La": "invmel",
    "Ce": "invmel",
    "Pb": "invmel",
    "Pr": "invmel",
    "Nd": "invmel",
    "Sr": "invmel",
    "Zr": "invmel",
    "Hf": "invmel",
    "Sm": "invmel",
    "Eu": "invmel",
    "Ti": "invmel",
    "Gd": "invmel",
    "Tb": "invmel",
    "Dy": "invmel",
    "Ho": "invmel",
    "Y": "invmel",
    "Er": "invmel",
    "Yb": "invmel",
    "Lu": "invmel",
}
"""
The method to use for each element when a method isn't otherwise specified.
"""

workman05_dmm = {
    "Rb": 0.05,
    "Ba": 0.563,
    "Th": 0.0079,
    "U": 0.0032,
    "Nb": 0.1485,
    "Ta": 0.0096,
    "La": 0.192,
    "Ce": 0.550,
    "Pb": 0.018,
    "Pr": 0.107,
    "Nd": 0.581,
    "Sr": 7.664,
    "Zr": 5.082,
    "Hf": 0.157,
    "Sm": 0.239,
    "Eu": 0.096,
    "Ti": 716.3,
    "Gd": 0.358,
    "Tb": 0.070,
    "Dy": 0.505,
    "Ho": 0.115,
    "Y": 3.328,
    "Er": 0.348,
    "Yb": 0.365,
    "Lu": 0.058,
}
"""
The trace element concentrations in the depleted MORB mantle from Workman & Hart (2005). All
concentrations are in ppmw.
"""

workman05_D = {
    "Rb": 1e-5,
    "Ba": 0.00012,
    "Th": 0.001,
    "U": 0.0011,
    "Nb": 0.0034,
    "Ta": 0.0034,
    "La": 0.01,
    "Ce": 0.022,
    "Pb": 0.014,
    "Pr": 0.027,
    "Nd": 0.031,
    "Sr": 0.025,
    "Zr": 0.033,
    "Hf": 0.035,
    "Sm": 0.045,
    "Eu": 0.050,
    "Ti": 0.058,
    "Gd": 0.056,
    "Tb": 0.068,
    "Dy": 0.079,
    "Ho": 0.084,
    "Y": 0.088,
    "Er": 0.097,
    "Yb": 0.115,
    "Lu": 0.120,
}
"""
The bulk partition coefficients for MORB production from Workman & Hart (2005).
"""

stracke03_bsic = {
    "Rb": 0.57,
    "Ba": 6.59,
    "Th": 0.088,
    "U": 0.027,
    "Nb": 1.95,
    "Ta": 0.124,
    "La": 1.68,
    "Ce": 5.89,
    "Pb": 0.09,
    "Nd": 7.45,
    "Sr": 81.0,
    "Zr": 64.0,
    "Hf": 1.78,
    "Sm": 2.69,
    "Eu": 1.04,
    "Ti": 7735.0,
    "Gd": 4.03,
    "Dy": 5.01,
    "Y": 28.5,
    "Er": 3.13,
    "Yb": 2.99,
    "Lu": 0.45,
}
"""
The trace element concentrations (ppmw) in bulk subducted igneous crust from Stracke et al. (2003).
"""

palme13_pm = {
    "Rb": 0.605,
    "Ba": 6.85,
    "Th": 0.0849,
    "U": 0.0229,
    "Nb": 0.595,
    "Ta": 0.043,
    "La": 0.6832,
    "Ce": 1.7529,
    "Pb": 0.185,
    "Pr": 0.2657,
    "Nd": 1.341,
    "Sr": 22.0,
    "Zr": 10.3,
    "Hf": 0.3014,
    "Sm": 0.4347,
    "Eu": 0.1665,
    "Ti": 1265.0,
    "Gd": 0.5855,
    "Tb": 0.1075,
    "Dy": 0.7239,
    "Ho": 0.1597,
    "Y": 4.13,
    "Er": 0.4684,
    "Yb": 0.4774,
    "Lu": 0.07083,
}
"""
The composition of the primitive mantle (ppmw) from Palme & O'Neill (2013).
"""

palme13_ci = {
    "Rb": 2.32,
    "Ba": 2.42,
    "Th": 0.03,
    "U": 0.00810,
    "Nb": 0.283,
    "Ta": 0.015,
    "La": 0.2414,
    "Ce": 0.6194,
    "Pb": 2.62,
    "Pr": 0.09390,
    "Nd": 0.4737,
    "Sr": 7.79,
    "Zr": 3.63,
    "Hf": 0.1065,
    "Sm": 0.1536,
    "Eu": 0.05883,
    "Ti": 447.0,
    "Gd": 0.2069,
    "Tb": 0.03797,
    "Dy": 0.2558,
    "Ho": 0.05644,
    "Y": 1.46,
    "Er": 0.1655,
    "Yb": 0.1687,
    "Lu": 0.02503,
}
"""
Trace element concentrations in a CI chondrite (ppmw) from Palme & O'Neill (2013).
"""

# From Gibson & Geist compilation
olv_D = {
    "Rb": 0.0003,
    "Ba": 0.000005,
    "Th": 0.00005,
    "U": 0.00038,
    "Nb": 0.0005,
    "Ta": 0.0005,
    "La": 0.0005,
    "Ce": 0.0005,
    "Pb": 0.003,
    "Pr": 0.0008,
    "Nd": 0.00042,
    "Sr": 0.00004,
    "Zr": 0.0033,
    "Hf": 0.0022,
    "Sm": 0.0011,
    "Eu": 0.0016,
    "Ti": 0.015,
    "Gd": 0.0011,
    "Tb": 0.0015,
    "Dy": 0.0027,
    "Ho": 0.0016,
    "Y": 0.0099,
    "Er": 0.013,
    "Yb": 0.020,
    "Lu": 0.020,
}
"""
Trace element partition coefficients between olivine and melt, compiled by Gibson & Geist (2010).
"""

# From Gibson & Geist compilation
opx_D = {
    "Rb": 0.0002,
    "Ba": 0.000006,
    "Th": 0.002,
    "U": 0.002,
    "Nb": 0.004,
    "Ta": 0.004,
    "La": 0.0031,
    "Ce": 0.0040,
    "Pb": 0.009,
    "Pr": 0.0048,
    "Nd": 0.01200,
    "Sr": 0.0007,
    "Zr": 0.013,
    "Hf": 0.03,
    "Sm": 0.0200,
    "Eu": 0.0130,
    "Ti": 0.086,
    "Gd": 0.0130,
    "Tb": 0.0190,
    "Dy": 0.0110,
    "Ho": 0.0065,
    "Y": 0.052,
    "Er": 0.045,
    "Yb": 0.080,
    "Lu": 0.120,
}
"""
Trace element partition coefficients between orthopyroxene and melt, compiled by Gibson & Geist
(2010).
"""

# From Gibson & Geist compilation
cpx_D = {
    "Rb": 0.0004,
    "Ba": 0.0004,
    "Th": 0.0059,
    "U": 0.0094,
    "Nb": 0.015,
    "Ta": 0.015,
    "La": 0.0490,
    "Ce": 0.0800,
    "Pb": 0.012,
    "Pr": 0.126,
    "Nd": 0.17800,
    "Sr": 0.091,
    "Zr": 0.119,
    "Hf": 0.284,
    "Sm": 0.2930,
    "Eu": 0.3350,
    "Ti": 0.350,
    "Gd": 0.3500,
    "Tb": 0.4030,
    "Dy": 0.4000,
    "Ho": 0.4270,
    "Y": 0.426,
    "Er": 0.420,
    "Yb": 0.400,
    "Lu": 0.376,
}
"""
Trace element partition coefficients between clinopyroxene and melt, compiled by Gibson & Geist
(2010).
"""

# From Gibson & Geist compilation
grt_D = {
    "Rb": 0.0002,
    "Ba": 0.00007,
    "Th": 0.009,
    "U": 0.028,
    "Nb": 0.015,
    "Ta": 0.015,
    "La": 0.0010,
    "Ce": 0.0050,
    "Pb": 0.005,
    "Pr": 0.014,
    "Nd": 0.05200,
    "Sr": 0.0007,
    "Zr": 0.270,
    "Hf": 0.400,
    "Sm": 0.2500,
    "Eu": 0.4960,
    "Ti": 0.600,
    "Gd": 0.84800,
    "Tb": 1.4770,
    "Dy": 2.2000,
    "Ho": 3.3150,
    "Y": 3.100,
    "Er": 4.400,
    "Yb": 6.600,
    "Lu": 7.100,
}
"""
Trace element partition coefficients between garnet and melt, compiled by Gibson & Geist (2010).
"""

# alphaMELTS defaults
spn_D = {
    "Rb": 0.0001,
    "Ba": 0.0001,
    "Th": 0.0,
    "U": 0.0,
    "Nb": 0.0,
    "Ta": 0.0,
    "La": 0.0100,
    "Ce": 0.0100,
    "Pb": 0.0,
    "Pr": 0.01,
    "Nd": 0.0100,
    "Sr": 0.0,
    "Zr": 0.0,
    "Hf": 0.0,
    "Sm": 0.0100,
    "Eu": 0.0100,
    "Ti": 0.15,
    "Gd": 0.0100,
    "Tb": 0.0100,
    "Dy": 0.0100,
    "Ho": 0.0100,
    "Y": 0.01,
    "Er": 0.0100,
    "Yb": 0.0100,
    "Lu": 0.0100,
}
"""
Trace element partition coefficients between spinel and melt, compiled by Gibson & Geist (2010).
"""

plg_D = {
    "Rb": 0.03,
    "Ba": 0.33,
    "Th": 0.05,
    "U": 0.11,
    "Nb": 0.01,
    "Ta": 0.0,
    "La": 0.2700,
    "Ce": 0.200,
    "Pb": 0.36,
    "Pr": 0.17,
    "Nd": 0.1400,
    "Sr": 2.0,
    "Zr": 0.01,
    "Hf": 0.01,
    "Sm": 0.1100,
    "Eu": 0.7300,
    "Ti": 0.04,
    "Gd": 0.0660,
    "Tb": 0.0600,
    "Dy": 0.0550,
    "Ho": 0.0480,
    "Y": 0.03,
    "Er": 0.0100,
    "Yb": 0.031,
    "Lu": 0.0250,
}
"""
Trace element partition coefficients between plagioclase and melt, compiled by Gibson & Geist
(2010).
"""

defaultD = _pd.DataFrame(
    {
        "olv": olv_D,
        "cpx": cpx_D,
        "opx": opx_D,
        "plg": plg_D,
        "grt": grt_D,
        "spn": spn_D,
    }
)
"""
Dataframe containing the partition coefficients for each phase in INVMEL.
"""

klb1_MineralProportions = _pd.DataFrame(
    [
        [0.609, 0.125, 0.119, 0.147, 0.000, 0.000],
        [0.597, 0.233, 0.158, 0.000, 0.012, 0.000],
        [0.646, 0.208, 0.076, 0.000, 0.000, 0.070],
    ],
    columns=["olv", "opx", "cpx", "grt", "spn", "plg"],
    index=["grt_field", "spn_field", "plg_field"],
)
"""
Mineral proportions (Wt%) for KLB1 in the garnet-, spinel-, and plagioclase-field (estimated
from Jennings and Holland, 2015).
"""

kg1_MineralProportions = _pd.DataFrame(
    [
        [0.181, 0.012, 0.422, 0.385, 0.000, 0.000],
        [0.110, 0.178, 0.641, 0.000, 0.071, 0.000],
        [0.118, 0.150, 0.655, 0.000, 0.000, 0.067],
    ],
    columns=["olv", "opx", "cpx", "grt", "spn", "plg"],
    index=["grt_field", "spn_field", "plg_field"],
)
"""
Mineral proportions (Wt%) for KG1 in the garnet-, spinel-, and plagioclase-field (estimated
from Matthews et al., 2021).
"""

mo91_MineralProportions = _pd.DataFrame(
    [
        [0.598, 0.221, 0.076, 0.115, 0.000, 0.000],
        [0.578, 0.270, 0.119, 0.000, 0.033, 0.000],
        [0.636, 0.263, 0.012, 0.000, 0.000, 0.089],
    ],
    columns=["olv", "opx", "cpx", "grt", "spn", "plg"],
    index=["grt_field", "spn_field", "plg_field"],
)
"""
Mineral proportions (Wt%) for lherzolite in the garnet-, spinel-, and plagioclase-field, from
McKenzie & O'Nions (1991).
"""


[docs]class species(object): """ The species base class. The class contains all the parameters and methods required to calculate the concentration of chemical species in melts, given T, P, F, etc. Parameters ---------- name : str The name of the species. c0 : float The concentration of the species in the solid. Must be the same units as required in the output. """ def __init__(self, name, c0, **kwargs): self.calculation_type = None self.name = name self.c0 = c0
[docs] def composition(self, state): """ Returns the concentration of the species in the melt for the specified state. This function should be redefined according to the chemical model being used. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. Returns ------- float The concentration of the species in the melt. """ return _np.nan
[docs]class batchSpecies(species): """ Implementation of the species class for batch melting with a constant partition coefficient. Parameters ---------- name : str The name of the species. c0 : float The concentration of the species in the solid. Must be the same units as required in the output. D : float The partition coefficient """ def __init__(self, name, c0, D, **kwargs): self.calculation_type = "accumulated" self.name = name self.c0 = c0 self._D = D
[docs] def composition(self, state): """ Returns the concentration in the melt during batch melting. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. """ return self.c0 / (self.D(state) * (1 - state["F"]) + state["F"])
[docs] def D(self, state): """ The partition coefficient. If a constant partition coefficient is used it will return that value. If a variable coefficient is used it will call the function to calculate it. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. Returns ------- float The partition coefficient """ if callable(self._D): d = self._D(state) else: d = self._D return d
[docs]class continuousSpecies_instantaneous(species): """ Implementation of the species class for batch melting with a constant partition coefficient. Parameters ---------- name : str The name of the species. c0 : float The concentration of the species in the solid. Must be the same units as required in the output. D : float The partition coefficient """ def __init__(self, name, c0, D, phi=0.005, **kwargs): self.calculation_type = "instantaneous" self.name = name self.c0 = c0 self._D = D self.phi = phi
[docs] def composition(self, state): """ Returns the instantaneous concentration in the melt during near-fractional (continuous) melting. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. """ D = self.D(state) exponent = (1 - self.phi) * (1 - D) / ((1 - self.phi) * D + self.phi) return self.c0 / ((1 - self.phi) * D + self.phi) * (1 - state.F) ** exponent
[docs] def D(self, state): """ The partition coefficient. If a constant partition coefficient is used it will return that value. If a variable coefficient is used it will call the function to calculate it. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. Returns ------- float The partition coefficient """ if callable(self._D): d = self._D(state) else: d = self._D return d
[docs]class continuousSpecies_accumulated(species): """ Implementation of the species class for batch melting with a constant partition coefficient. Parameters ---------- name : str The name of the species. c0 : float The concentration of the species in the solid. Must be the same units as required in the output. D : float The partition coefficient """ def __init__(self, name, c0, D, phi=0.005, **kwargs): self.calculation_type = "accumulated" self.name = name self.c0 = c0 self._D = D self.phi = phi
[docs] def composition(self, state): """ Returns the concentration in the melt during batch melting. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. """ Deff = (1 - self.phi) * self.D(state) + self.phi return self.c0 / state.F * (1 - (1 - state.F) ** (1 / Deff))
[docs] def D(self, state): """ The partition coefficient. If a constant partition coefficient is used it will return that value. If a variable coefficient is used it will call the function to calculate it. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. Returns ------- float The partition coefficient """ if callable(self._D): d = self._D(state) else: d = self._D return d
[docs]class mineralTransition(object): """ Template for the mineralTransition class used to define mineral transitions in the invMEL chemistry model (and perhaps other models in the future). Parameters ---------- parameters : dict A dictionary of parameters required by the transition function. """ def __init__(self, parameters={}): self.parameters = parameters
[docs] def transition(self, T): """ Template transition method. Must be a function of temperature using parameters contained in the parameters dictionary. Parameters ---------- T : float Temperature (in degC). Returns ------- float Transition pressure (in GPa). """ return _np.nan
[docs]class mineralTransition_linear(mineralTransition): """ Representation of a mineral transition (e.g. garnet-in) as a linear function of temperature. The parameters dictionary must contain items: - 'intercept': the intercept pressure (T=0) in GPa - 'gradient': the gradient of the line (GPa/degC) """
[docs] def transition(self, T): return self.parameters["intercept"] + self.parameters["gradient"] * T
[docs]class mineralTransition_isobaric(mineralTransition): """ Representation of a mineral transition that does not depend on temperature. Temperature is still taken as an argument for uniformity. The parameters dictionary must contain items: - 'transition_pressure': the transition pressure in GPa """
[docs] def transition(self, T): return self.parameters["transition_pressure"]
[docs]class invmelSpecies(species): """ Implementation of the forward trace element model used by invmel (McKenzie & O'Nions, 1991). The default values of the mineral-melt partition coefficients are taken from REFERENCE. Parameters ---------- name : str The name of the species. c0 : float The concentration of the species in the solid. Must be the same units as required in the output. olv_D : float The partition coefficient for olivine-melt. cpx_D : float The partition coefficient for clinopyroxene-melt. opx_D : float The partition coefficient for orthopyroxene-melt. spn_D : float The partition coefficient for spinel-melt. grt_D : float The partition coefficient for garnet-melt. plg_D : float The partition coefficient for plagioclase-melt. MineralProportions : pandas.DataFrame, default: pyMelt.chemistry.mo91_MineralProportions A dataframe with the proportions of each mineral phase (columns) in the garnet-, spinel-, and plagioclase-field for the lithology. See built in defaults for formatting of keys and columns: mo91_MineralProportions, klb1_MineralProportions, kg1_MineralProportions. density : float, default: 3.3 The density of the mantle (g cm-3) cpxExhaustion : int, default: 0.18 The melt fraction at which cpx (and grt/plg/spn) are exhausted. garnetOut : mineralTransition, default: 1/666.7*T + 400/666.7 GPa The garnet-out reaction (decompression). spinelIn : mineralTransition, default: 1/666.7*T + 533/666.7 GPa The spinel-in reaction (decompression). spinelOut : mineralTransition, default: 25*0.033 GPa The spinel-out reaction (decompression). plagioclaseIn : mineralTransition, default: 35*0.033 GPa The plagioclase-in reaction (decompression). Attributes ---------- name : str Name of the species c0 : float Concentration of the species in the solid. D : numpy.Array The partition coefficients in the order: olivine, cpx, opx, spinel, garnet, plag. MineralProportions_solid : pandas.DataFrame See parameters above. density : float The density of the mantle (g cm-3) """ def __init__( self, name, c0, olv_D, cpx_D, opx_D, spn_D, grt_D, plg_D, mineralProportions=mo91_MineralProportions, density=3.3, cpxExhaustion=0.18, garnetOut=mineralTransition_linear({'gradient': 1 / 666.7, 'intercept': 400 / 666.7}), spinelIn=mineralTransition_linear({'gradient': 1 / 666.7, 'intercept': 533 / 666.7}), spinelOut=mineralTransition_isobaric({'transition_pressure': 25.0 * 0.033}), plagioclaseIn=mineralTransition_isobaric({'transition_pressure': 35.0 * 0.033}), **kwargs ): self.calculation_type = "instantaneous" self.name = name self.c0 = c0 self._D = { "olv": olv_D, "cpx": cpx_D, "opx": opx_D, "spn": spn_D, "grt": grt_D, "plg": plg_D, } self.mineralProportions_solid = mineralProportions self.density = density self.cpxExhaustion = cpxExhaustion self.garnetOut = garnetOut self.spinelIn = spinelIn self.spinelOut = spinelOut self.plagioclaseIn = plagioclaseIn self._cs = c0 self._F_prev = 0.0 self._cl_prev = None
[docs] def composition(self, state): # Check if this is a new calculation or not: if state.F < self._F_prev: self._cs = self.c0 if state.F == 1: # If the lithology is immediately fully molten: if self._cl_prev is None: return self._cs else: return self._cl_prev D = self.D_bulk(state["P"], state["T"], state.F) Pbar = self.P_bulk(state["P"], state["T"], state.F) if D < 1e-4: _warn( self.name + " is extremely incompatible, unless the step size is extremely small" " its partitioning behaviour is unlikely to be captured correctly. " "You are probably better off calculating it using the " "ContinuousSpecies_accumulated class." ) k1 = self._dcsdX(self._F_prev, self._cs, self._cl(self._cs, state.F, D, Pbar), D, Pbar) k2 = self._dcsdX( self._F_prev + (state.F - self._F_prev) / 2, self._cs + k1 * (state.F - self._F_prev) / 2, self._cl(self._cs + k1 * (state.F - self._F_prev) / 2, self._F_prev + (state.F - self._F_prev) / 2, D, Pbar), D, Pbar, ) k3 = self._dcsdX( self._F_prev + (state.F - self._F_prev) / 2, self._cs + k2 * (state.F - self._F_prev) / 2, self._cl(self._cs + k2 * (state.F - self._F_prev) / 2, self._F_prev + (state.F - self._F_prev) / 2, D, Pbar), D, Pbar, ) k4 = self._dcsdX( self._F_prev + (state.F - self._F_prev), self._cs + k3 * (state.F - self._F_prev), self._cl(self._cs + k3 * (state.F - self._F_prev), self._F_prev + (state.F - self._F_prev), D, Pbar), D, Pbar, ) cs = self._cs + (1 / 6) * (state.F - self._F_prev) * (k1 + 2 * k2 + 2 * k3 + k4) cl = self._cl(cs, state.F, D, Pbar) # Check if discretisation is too coarse if (k1 + 2 * k2 + 2 * k3 + k4) > 0 and D < 1: _warn( "Discretisation is too coarse to capture the behaviour of " + self.name + "." ) cl = _np.nan self._cs = _np.nan # Prevent float errors elif cs < 1e-6: self._cs = 0 else: self._cs = cs self._F_prev = state.F self._cl_prev = cl return cl
[docs] def D(self, state): """ The partition coefficient. If a constant partition coefficient is used it will return that value. If a variable coefficient is used it will call the function to calculate it. Parameters ---------- state : pandas.Series The state of the system, e.g. temperature (T), pressure (P), melt fraction (F). This will most likely be generated automatically by the `MeltingColumn_1D` class. Returns ------- float The partition coefficient """ d = {} for min in self._D: if callable(self._D[min]): d[min] = self._D[min](state) else: d[min] = self._D[min] return d
def mineralProportions(self, P, T): GarnetSpinel = self._GarnetSpinelTransition(P, T) SpinelPlagioclase = self._SpinelPlagioclaseTransition(P, T) if GarnetSpinel == 1 and SpinelPlagioclase == 1: mineralProportions = self.mineralProportions_solid.loc["grt_field"] elif GarnetSpinel == 0 and SpinelPlagioclase == 1: mineralProportions = self.mineralProportions_solid.loc["spn_field"] elif GarnetSpinel == 0 and SpinelPlagioclase == 0: mineralProportions = self.mineralProportions_solid.loc["plg_field"] elif SpinelPlagioclase == 1: grtField = self.mineralProportions_solid.loc["grt_field"] spnField = self.mineralProportions_solid.loc["spn_field"] mineralProportions = [ (grtFieldProp * GarnetSpinel + spnFieldProp * (1 - GarnetSpinel)) for grtFieldProp, spnFieldProp in zip(grtField, spnField) ] mineralProportions = _pd.Series( mineralProportions, index=self.mineralProportions_solid.columns ) else: spnField = self.mineralProportions_solid.loc["spn_field"] plgField = self.mineralProportions_solid.loc["plg_field"] mineralProportions = [ (spnFldPrp * SpinelPlagioclase + plgFldPrp * (1 - SpinelPlagioclase)) for spnFldPrp, plgFldPrp in zip(spnField, plgField) ] mineralProportions = _pd.Series( mineralProportions, index=self.mineralProportions_solid.columns ) return mineralProportions def D_bulk(self, P, T, F=0.0): mineralProportions = self.mineralProportions(P, T) Dminerals = self.D(_pd.Series({"P": P, "T": T})) cpxExhaustion = self.cpxExhaustion if F < cpxExhaustion: D = sum([Dminerals[min] * mineralProportions[min] for min in mineralProportions.index]) else: D = ((Dminerals["olv"] * mineralProportions["olv"] + Dminerals["opx"] * mineralProportions["opx"]) / (mineralProportions["olv"] + mineralProportions["opx"])) return D def P_bulk(self, P, T, F=0.0): mineralProportions = self.mineralProportions(P, T) Dminerals = self.D(_pd.Series({"P": P, "T": T})) cpxExhaustion = self.cpxExhaustion p = {} if F < cpxExhaustion: p["cpx"] = mineralProportions["cpx"] / cpxExhaustion p["grt"] = mineralProportions["grt"] / cpxExhaustion p["spn"] = mineralProportions["spn"] / cpxExhaustion p["plg"] = mineralProportions["plg"] / cpxExhaustion p["olv"] = ( mineralProportions["olv"] * (1 - (p["cpx"] + p["grt"] + p["spn"] + p["plg"])) / (mineralProportions["olv"] + mineralProportions["opx"]) ) p["opx"] = ( mineralProportions["opx"] * (1 - (p["cpx"] + p["grt"] + p["spn"] + p["plg"])) / (mineralProportions["olv"] + mineralProportions["opx"]) ) else: p["olv"] = mineralProportions["olv"] / ( mineralProportions["olv"] + mineralProportions["opx"] ) p["opx"] = mineralProportions["opx"] / ( mineralProportions["olv"] + mineralProportions["opx"] ) p["cpx"] = 0 p["grt"] = 0 p["spn"] = 0 p["plg"] = 0 Pbar = sum([Dminerals[min] * p[min] for min in Dminerals.keys()]) return Pbar def _SpinelPlagioclaseTransition(self, P, T): """ Calculate proportions contributing to sp pl transition as contribution from the sp field mantle. Parameters ---------- P : float pressure in GPa T : float temperature in degC Returns ------- float contribution from sp mantle """ plagIn = self.plagioclaseIn.transition(T) spinelOut = self.spinelOut.transition(T) if P <= spinelOut: SpinelPlagioclaseTransition = 0 elif P >= plagIn: SpinelPlagioclaseTransition = 1 else: SpinelPlagioclaseTransition = 1 - (plagIn - P) / (plagIn - spinelOut) return SpinelPlagioclaseTransition def _GarnetSpinelTransition(self, P, T): """ Calculate proportions contributing to gt sp transition as contribution from the gt field mantle. Parameters ---------- P : float pressure in GPa T : float temperature in degC Returns ------- float contribution from gt mantle """ GarnetOut = self.garnetOut.transition(T) SpinelIn = self.spinelIn.transition(T) if P >= SpinelIn: GarnetSpinelTransition = 1 elif P <= GarnetOut: GarnetSpinelTransition = 0 else: GarnetSpinelTransition = 1 - (P - SpinelIn) / (GarnetOut - SpinelIn) return GarnetSpinelTransition def _dcsdX(self, X, cs, cl, Dbar, Pbar): """ Rate of change of rare-earth element concentration in point average solid. Parameters ---------- X : float melt fraction cs : float point average solid residue composition cl : float liquid composition Dbar : float bulk distribution coefficient for solid assemblage Pbar : float bulk distribution coefficient for melting package Returns ------- float rate of change of REE concentration in point average solid residue with respect to melt fraction """ dcsdX = (cs - cl) / (1 - X) return dcsdX def _cl(self, cs, X, Dbar, Pbar): """ Calculates instantaneous melt composition generated from a point average solid. Parameters ---------- X : float melt fraction cs : float point average solid residue composition Dbar : float bulk distribution coefficient for solid assemblage Pbar : float bulk distribution coefficient for melting package Returns ------- float instantaneous melt composition """ cl = cs * (1 - X) / (Dbar - Pbar * X) return cl