MagmaForge: Equilibrium Crystallization by heat extraction

When you crystallize a near-eutectic liquid, much of the crystallization occurs over a very narrow temperature range. Driving the calculation by extracting energy directly, rather than by fixing the temperature provides an alternative way of exploring this phenomenon.

Note that this is a prototype notebook and may show some anomalous behaviour for some combinations of bulk composition and pressure.

Initialization

Import necessary packages:

from thermoengine import magmaforge

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Define a bulk composition:

# Early Bishop Tuff B7-1 (Anderson et al., 2000)
comp_wtpt = pd.Series({
    'SiO2': 77.6,
    'Al2O3': 12.6,
    'FeO': 0.67,
    'MgO': 0.02,
    'CaO': 0.44,
    'Na2O': 4.15,
    'K2O': 4.67,
    'H2O': 7.0
})

System and Calculations

Set the pressure and entropy step for the calculation:

P_bar = 1000.0
# P_bar = 1800.0 # This is very close to the three-phase saturation pressure. It will take a long time to run!
dS = 0.25 # J/K

Run the calculation:

sys = magmaforge.System(comp=comp_wtpt,
                        T_liquidus=True,
                        # T_C=825.0,
                        P_bar=P_bar,
                        Fe3_tFe = 0.8,
                        # logfO2=('QFM', 1.0),
                        database='MELTS_v1_0'
                        )

S_start = sys.total_entropy

sys.crystallize_by_entropy_reduction(S_step=dS)

<thermoengine.magmaforge.system.System object at 0xffff8f1660b0>

Calculation Results

 Make a plot of phase abundances vs heat extracted and temperature (for comparison):

# Assemble calculation results:
phase_fracs = sys.history.get_phase_frac_table(phases='present')
entropy_steps = sys.history.get_total_entropies()
heat_extracted = (S_start - entropy_steps)*sys.history.get_temperatures()/1000.0

fig, ax = plt.subplots(1,2, sharey='row')

line_to_plot = np.zeros(len(phase_fracs))
line_to_plot_prev = line_to_plot.copy()
phases_present = list(phase_fracs.columns)
phases_to_plot = []
for phs in phases_present:
    if phs != 'Liquid' and phs not in phases_to_plot:
        phases_to_plot.append(phs)
        if phs+'_1' in phases_present:
            phases_to_plot.append(phs+'_1')

for phase in phases_to_plot:
    line_to_plot += phase_fracs[phase]
    ax[0].fill_between(heat_extracted,
                       line_to_plot*100,
                       line_to_plot_prev*100,
                       label=phase)
    ax[1].fill_between(phase_fracs.index - 273.15,
                       line_to_plot*100,
                       line_to_plot_prev*100,
                       label=phase)

    line_to_plot_prev = line_to_plot.copy()

ax[0].legend(fontsize=8)
ax[0].set_xlabel('Heat Extracted (kJ)')
ax[0].set_xlim(np.max(heat_extracted), np.min(heat_extracted))
ax[1].set_xlim(phase_fracs.index[-1] - 273.15, phase_fracs.index[0] - 273.15)
ax[1].set_xlabel('Temperature (˚C)')
for a in ax:
    a.set_ylim(0,100)
ax[0].set_ylabel('Phase Fraction (%)')

fig.tight_layout()

plt.show()

Plot temperature vs heat extracted:

fig, ax = plt.subplots()

ax.plot(heat_extracted, phase_fracs.index - 273.15)
ax.set_xlabel('Heat Extracted (kJ)')
ax.set_ylabel('Temperature (˚C)')

Text(38.347222222222214, 0.5, 'Temperature (˚C)')