-
Notifications
You must be signed in to change notification settings - Fork 41
/
partitioning.py
64 lines (49 loc) · 2.23 KB
/
partitioning.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
# the Earth and Planetary Sciences
# Copyright (C) 2012 - 2019 by the BurnMan team, released under the GNU
# GPL v2 or later.
from __future__ import absolute_import
import numpy as np
from .. import constants
def calculate_nakajima_fp_pv_partition_coefficient(
pressure, temperature, bulk_composition_mol, initial_distribution_coefficient
):
"""
Calculate the partitioning of iron between periclase and bridgmanite as given
by Nakajima et al., 2012.
:param pressure: Equilibrium pressure [Pa]
:type pressure: float
:param temperature: Equilibrium temperature [K]
:type temperature: float
:param bulk_composition_mol: Bulk composition [mol].
Only Mg, Fe, and Si are assumed to explicitly affect the partitioning
with Al playing an implicit role.
:type bulk_composition_mol: dict
:param initial_distribution_coefficient: The distribution coefficient (Kd_0)
at 25 GPa and 0 K.
:type initial_distribution_coefficient: float
:returns: The proportion of Fe in ferropericlase and perovskite, respectively
:rtype: tuple
"""
norm = bulk_composition_mol["Mg"] + bulk_composition_mol["Fe"]
f_FeO = bulk_composition_mol["Fe"] / norm
f_SiO2 = bulk_composition_mol["Si"] / norm
Kd_0 = initial_distribution_coefficient
delV = 2.0e-7 # in m^3/mol, average taken from Nakajima et al 2012, JGR
# eq 5 Nakajima et al 2012, JGR. Solved for ln(K(P,T,X))
rs = ((25.0e9 - pressure) * delV / (constants.gas_constant * temperature)) + np.log(
Kd_0
)
# The exchange coefficent at P and T. K(P,T,X) in eq 5 Nakajima et al 2012
K = np.exp(rs)
# Solving equation 6 in Nakajima et al., 2012 for X_Fe_fp and X_Fe_pv
# Solved using the definition of the distribution coefficient
# to define X_Fe_fp as a function of X_Fe_pv
num_to_sqrt = (-4.0 * f_FeO * (K - 1.0) * K * f_SiO2) + (
pow(1.0 + (f_FeO * (K - 1)) + ((K - 1.0) * f_SiO2), 2.0)
)
X_Fe_pv = (
-1.0 + f_FeO - (f_FeO * K) + f_SiO2 - (f_SiO2 * K) + np.sqrt(num_to_sqrt)
) / (2.0 * f_SiO2 * (1.0 - K))
X_Fe_fp = X_Fe_pv / (((1.0 - X_Fe_pv) * K) + X_Fe_pv)
return (X_Fe_fp, X_Fe_pv)