MagmaSat: Plot Isobars and Isopleths

Calculates the CO₂-H₂O concentrations in melt in equilibrium with a CO₂-H₂O mixed fluid at specified pressure and temperature.

This uses the VESIcal python library and is a preview of a forthcoming version of VESIcal compatible with the new version of thermoengine.

Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):

To install the preview version of VESIcal that works with the new pip installable version of thermoengine, run this command (uncomment first):

%pip install kaylai/VESIcal@thermoenginelite

import VESIcal as v

Specify the composition of the liquid. Here a basaltic composition is given as an example:

oxide_comp = {
        'SiO2':49.91,
        'TiO2': 1.47,
        'Al2O3': 17.91,
        'Fe2O3': 2.45,
        'FeO':  7.02,
        'MnO':  0.16,
        'MgO':  6.62,
        'CaO':  10.02,
        'Na2O':  3.02,
        'K2O':  0.64,
        'P2O5':  0.2,
        'H2O': 0.1,
        'CO2': 0.1
        }

Now convert this to a VESIcal Sample object (this allows the code to deal with unit conversions, normalisation, etc)

sample = v.Sample(oxide_comp)

Do the calculation

VESIcal will iterate over fluid compositions (from pure H$_2$O to pure CO$_2$) and pressure to identify how much H$_2$O and CO$_2$ will be dissolved in the melt. The results are isobars (constant pressure, varying fluid composition) and isopleths (constant fluid composition, varying pressure).

isobars, isopleths = v.calculate_isobars_and_isopleths(sample=sample,
                                            temperature=1000.0, # degC
                                            pressure_list=[1000.0, 2000.0, 3000.0], # bars
                                            isopleth_list=[0.25,0.5,0.75], # Mole fraction
                                            print_status=True,
                                            ).result

Calculating isobar at 1000.0 bars

 Calculating isobar control point at XH2Ofluid = 0
 Calculating isopleth at XH2Ofluid = 0.25
 Calculating isopleth at XH2Ofluid = 0.5
 Calculating isopleth at XH2Ofluid = 0.75
 Calculating isobar control point at XH2Ofluid = 1
 done.
Calculating isobar at 2000.0 bars

 Calculating isobar control point at XH2Ofluid = 0
 Calculating isopleth at XH2Ofluid = 0.25
 Calculating isopleth at XH2Ofluid = 0.5
 Calculating isopleth at XH2Ofluid = 0.75
 Calculating isobar control point at XH2Ofluid = 1
 done.
Calculating isobar at 3000.0 bars

 Calculating isobar control point at XH2Ofluid = 0
 Calculating isopleth at XH2Ofluid = 0.25
 Calculating isopleth at XH2Ofluid = 0.5
 Calculating isopleth at XH2Ofluid = 0.75
 Calculating isobar control point at XH2Ofluid = 1
 done.
Done!

Make a plot!

 VESIcal provides a built-in method for visualising these results quickly:

v.plot(isobars=isobars, isopleths=isopleths)

(<Figure size 1200x800 with 1 Axes>, <Axes: xlabel='H$_2$O wt%', ylabel='CO$_2$ wt%'>)

but you can also make a custom plot by interacting with the isobars and isopleths variables directly. You could save them to a csv file, e.g.:

isobars.to_csv('my_isobars.csv')

# or plot them directly using a plotting library like `matplotlib`.