PEC model

Start by importing MagmaPEC and MagmaPandas and any other packages you want to use. Here we also import Pandas for importing pressure data. For details on the use of MagmaPandas, please see it’s documentation.

import MagmaPEC as mpc
import MagmaPandas as mp

import pandas as pd

Import your melt inclusion and olivine data.

melt_file = "./data/melt.csv"
olivine_file = "./data/olivine.csv"

melt = mp.read_melt(melt_file, index_col=["name"])
olivine = mp.read_olivine(olivine_file, index_col=["name"])

If you have data for internal inclusion pressures (in bars) you can also import those.

pressure_file ="./data/pressure.csv"

pressure = pd.read_csv(pressure_file, index_col = ["name"]).squeeze()

pressure

name
PI032-04-01    6527.087712
PI032-04-02    6828.021536
PI041-02-02    4802.872035
PI041-03-01    6958.250380
PI041-03-03    6826.596778
PI041-05-04    4013.630505
PI041-05-06    5269.481653
PI041-07-01    4969.100477
PI041-07-02    4175.590768
PI052-01-02    2388.318094
Name: 0, dtype: float64

Alternatively, if you have melt CO2 (and H2O) data, you can use MagmaPandas to calculate volatile saturation pressures. First calculate temperature with some initial pressure and then iteratively calculate pressure and temperature. See the MagmaPandas documentation for more details.

temperature = melt.temperature(P_bar = 5e3)
pressure = melt.volatile_saturation_pressure(T_K=temperature)

while True:
    temperature = melt.temperature(P_bar=pressure)
    pressure_new = melt.volatile_saturation_pressure(T_K=temperature)
    dP = (pressure_new - pressure) / pressure_new # calculate percentage change
    pressure = pressure_new.copy()
    if (dP < 0.01).all():
        # break when all pressure have converged within 1% of previous values.
        break

pressure

Saturation pressure...    |█████████████████████████| 100% [10/10] in 0.9s
Saturation pressure...    |█████████████████████████| 100% [10/10] in 0.9s

name
PI032-04-01    6156.914986
PI032-04-02    6828.021538
PI041-02-02    4802.872035
PI041-03-01    6958.250380
PI041-03-03    6826.596776
PI041-05-04    4013.630505
PI041-05-06    5269.481651
PI041-07-01    4969.100476
PI041-07-02    4175.590771
PI052-01-02    2388.318094
dtype: float64

Here we will also use a whole-rock dataset to set up a model for predicting the initial FeO contents of our melt inclusions. For explanation of this model, please follow the Initial FeO prediction example.

wholerock_file = "./data/wholerock.csv"

wholerock = mp.read_melt(wholerock_file, index_col=["name"])

Let’s begin by setting up the model to predict the initial FeO content of our melt inclusions according to the previous example. Here, the FeO prediction model is based on TiO2, Al2O3 and CaO, which we can check by accessing their coefficients:

x = wholerock.drop(columns=["FeO"])
FeOi_predict = mpc.FeOi_prediction(x=x, FeO=wholerock["FeO"])

do_not_use = ["MnO", "P2O5", "Cr2O3", "total"]

model_fits = FeOi_predict.calculate_model_fits(exclude=do_not_use)
FeOi_predict.select_predictors(idx=3)

FeOi_predict.coefficients

To use this model during PEC correction, we need to store this model as a callable function, which we can do with the model attribute:

FeO_model = FeOi_predict.model

Let’s first quickly test the model by predicting FeO contents for our uncorrected melt compositions:

FeO_model(melt)

This all looks good, so let continue with setting up the PEC correction model.

First, preview the melt and olivine data to check everything looks ok:

melt.head()

olivine.head()

Make sure that each row in melt and olivine is a matching pair of melt inclusion and olivine host. Here we use the sample names stored in the indices of both dataframes:

print(melt.index, "\n\n", olivine.index)

We can verify that the indices are the same with the equals method:

melt.index.equals(olivine.index)

Check the configuration of MagmaPandas and the PEC model and make changes if you want to (see the configuration example). You can change MagmaPandas settings either via the MagmaPandas object (here: mp.configuration) or the MagmaPEC object (here: mpc.configuration).

print(mpc.configuration)
print(mpc.PEC_configuration)

Now we can initialise the PEC model, where we need the following data:

• inclusions

  inclusion major element compositions in oxide wt. % as a MagmaPandas Melt frame.

• olivines

  olivine major element compositions in oxide wt. % as a MagmaPandas Olivine frame.

• P_bar

  pressures in bar at which to run the model. You can use a fixed pressure for all inclusions, e.g. P_bar=2e3 for 2 kbar, or indicate specific pressures per inclusion. In this example we do the latter, by passing the above imported pandas Series with internal inclusion pressures

• FeO_target

  Estimated initial FeO contents of melt inclusions. You can use a fixed value for all melt inclusions, e.g. FeO = 11, specific values for individual melt inclusions, stored in a pandas Series, or a predictive equation based on melt major element composition. This equation needs to be a callable function that accepts a Pandas DataFrame with melt compositions in oxide wt. % as input and return an array-like object with initial FeO contents per inclusion. In the example above we showed how to set up a function like this using MagmaPEC.

pec_model = mpc.PEC(inclusions=melt, olivines=olivine, P_bar=pressure, FeO_target=FeO_model)

Now we simply run the model with the correct method. It runs in two stages: the first stage equilibrates Fe and Mg between melt inclusions and olivine hosts through Fe-Mg cation exchange and the second stage corrects melt inclusions back to their initial FeO contents by melting or crystallising olivine.

This method returns three objects:

• corrected melt compositions as a MagmaPandas Melt frame

• PEC extents as a pandas DataFrame

• inclusion entrapment temperatures as a pandas Series

melts_corrected, pec, temperatures = pec_model.correct()

The dataframe with pec results has three columns:

• equilibration_crystallisation:

  crystallisation extents during the equilibration stage

• PE_crystallisation:

  crystallisation extents during the crystallisation stage

• total_crystallisation:

  total amount of post-entrapment crystallisation

with all data in percent.

pec

Results are also stored in the pec object and can be accessed via the inclusions and olivine_corrected attributes, while entrapment temperatures can be calculated on the fly with the temperature method of the inclusions Melt frame (see the MagmaPandas documentation)

pec_model.inclusions

pec_model.olivine_corrected

pec_model.inclusions.temperature(P_bar=pressure)