Problem with a phase equilibrium of a mineral assemblage

BurnMan
,
1

Dear community,
i have a problem in creating the correct python code to create the equilibrium of a mineral assemblage. I would like to reproduce the P-T conditions of this chemical equation (mineral assemblage): Dolomite + Enstatite = Diopside + Olivine + Garnet + CO2

Mi questions are:

  • how can i create the mineral assemblage?
  • how can i find the phase equilibrium of the equation at determined (experimentally) P-T conditions? And finally compare the quantum mechanical results with experimental ones?
  • can i use also the CO2 as a fluid in the equation?

Thank in advance to all. If you can help me i really thank you

2

Dear Carnivale,

We can certainly help you!

  • Equilibrium problems are described in detail in the 5th part of the BurnMan Tutorial: Part 5: Equilibrium problems - BurnMan 1.3.0a0 documentation
    There you can find examples of how to build an assemblage from different minerals, and select appropriate constraints (in your case, pressure and temperature).
  • Comparing the thermodynamic model results (not quantum mechanical) with experimental ones is up to you. Which things do you want to compare?
  • A CO2 fluid object can be initialised from CO2 = burnman.minerals.HP_2011_fluids.CO2()

It isn’t clear from your question exactly what you want to do.

  • Do you want to find the position of the P-T reaction line for an endmember reaction? e.g. dol + di = fo + cc + CO2? Your equation is not balanced (e.g. there’s no Al2O3 on the LHS)
  • Do you want to find the equilibrium state for an assemblage of known composition containing carbonate + orthopyroxene + clinopyroxene + olivine + garnet + fluid? If so, there are several independent reactions which will all take place in the rock at the same time, and it isn’t necessary to write any of them out (although BurnMan can also calculate all the potential balanced reactions for you!).

The tutorial is designed to show you how to solve both kinds of problems.

3

P.S. If you need any more help, please attach a file with what you’ve tried and explain where you are stuck.

Best wishes,
Bob

4

Hello Bob, first of all thank you very much for your kindly reply.

My problem is apparently easy to solve.
The equation mentioned before (Dolomite+Enstatite=Diopside+Olivine+Garnet+CO2) is a reaction that occurs under mantle conditions and there is an univariant point experimentally found where all mineral phases are in equilibrium (attached the figure from bibliography, Dalton and Presnall, 1998).

My idea is try to reproduce the experimentally conditions (P=2.6GPa, T=1260°C) and the same mineral assemblage, and see what P-T conditions i can find with those conditions using Burnman.

I know how to create the assemblage.
But then i have problem to find the Equilibrium Temperature for the entire assemblage.
See for example with fixed assemblage compositions and Pressure how the equilibrium temperature vary and compare the results with experimentally ones from literature.

The final goal is to depict a P-T diagram of the minerals assemblage with the Burnman results (such as in the first Figure example of the manual chapter 3.5.3).
Figure from Dalton and Presnall.pdf (108.4 KB)

5

This is the script that i used:

import numpy as np
import matplotlib.pyplot as plt
import burnman
from burnman import equilibrate
from burnman import minerals
Fo = minerals.SLB_2011.forsterite()
Di = minerals.SLB_2011.diopside()
En = minerals.SLB_2011.enstatite()
Gt = minerals.HHPH_2013.py()
Dol = minerals.HP_2011_ds62.dol()
Mst = minerals.HP_2011_ds62.mst()
CO2 = minerals.HP_2011_fluids.CO2()
from burnman.tools.chemistry import equilibrium_temperature
assemblage = burnman.Composite(phases=[Fo, Di, En, Gt, Dol, Mst, CO2],fractions=[0.2, 0.2, 0.2, 0.1, 0.1, 0.1, 0.1])
P = 2.6e9
T = equilibrium_temperature([En, Dol, Fo, Di, Gt, CO2], [78.0, 22.0, -36.0, -49.0, -5.0, -10.0], P, 1500)

Doesn’t work

The mineral assemblage (left), stoichiometry, pressure, temperature initial guess (right)

6

Ok, a few points:

  • Don’t mix and match endmembers from different datasets. SLB_2011 and HP_2011_ds62 are not compatible with each other. Use fo, di, en, py from the HP datasets.
  • You haven’t yet answered my question: do you want to look at an endmember reaction or an equilibrium involving solutions? The Dalton and Presnall paper includes both.
  • You have attempted to solve an endmember reaction problem, but your stiochiometries are not consistent. There is no Al2O3 in enstatite and dolomite, but there is in pyrope.
  • Your script suggests you are trying to recreate reaction 1A from the paper. This is not an endmember reaction, although it does look like one. opx, cpx and garnet in this equation are all solid solutions.

Endmember problems are solved using the strategy you tried; you just need to use the correct endmembers and the correct reaction. Here are the required tweaks:

fo = minerals.HP_2011_ds62.fo()
di = minerals.HP_2011_ds62.di()
en = minerals.HP_2011_ds62.en()
dol = minerals.HP_2011_ds62.dol()
CO2 = minerals.HP_2011_fluids.CO2()
assemblage = burnman.Composite([fo, di, en, dol, CO2])
reaction_stoichiometry = reactions_from_formulae(assemblage.endmember_formulae,
assemblage.endmember_names,
return_strings=False)
P = 2.6e9
T = equilibrium_temperature(assemblage.phases, reaction_stoichiometry[0], P, 1500)

The problem including solutions would also be easy to solve in burnman, but your problem just happens to need something that I haven’t implemented yet - the entropy of CORK EoS fluids. If you need this, let me know. You’ll have to give me a week or so to find time to write the necessary functions.

Best wishes,
Bob

7

P.S. dolomite and olivine are also solutions in that paper (although olivine is almost pure forsterite). I haven’t coded these up in burnman yet, so we would need to do that as well.

8

The final goal is try to look at both endmember reaction and equilibrium solutions.
It would be perfect if we include equilibrium involving solid solutions.

I don’t need it shortly, even better wait some weeks (I will be away on business until the beginning of the new year). So if you are available we can solve the problem on January.

We can avoid the Garnet and consider the reaction: CaMg(CO3)2 (dolomite) + 4MgSiO3 (Enstatite) = CaMgSi2O6 (Diopside) + 2Mg2SiO4 (Olivine) + 2CO2 (1)

The initial idea was to compute the thermodynamic of the equation (1) with Burnman and compare the results with the previous ones experimentally found.

Again thank you very much for you professional treatment. If you are interested we can also start a discussion and why not a consideration in our work.

Regards,
Gabriele

9

Hi Gabriele,

Ok, great, thanks for the information!

The additions to the code shouldn’t take long, I just need to find a bit of time. In checking your problem, I also spotted another new feature that would make coding up your problem easier. You can track and comment on the issues here: Preserve phase composition when transforming solution basis · Issue #559 · geodynamics/burnman · GitHub and Implement derivative properties in CORK EoS · Issue #560 · geodynamics/burnman · GitHub.

If you wanted to help with code development, I would be very happy to help with code review / advice. It’s not necessary (I’m happy to make the tweaks), but if you contributed code and documentation I’d be happy to add you to the contributor/author list.

We can avoid the Garnet and consider the reaction:
CaMg(CO3)2 (dolomite) + 4MgSiO3 (Enstatite) = CaMgSi2O6 (Diopside) + 2Mg2SiO4 (Olivine) + 2CO2 (1)

Excellent, this is the equation I used above. The line:

reaction_stoichiometry = reactions_from_formulae(assemblage.endmember_formulae,
assemblage.endmember_names,
return_strings=False)

automatically calculated all possible reaction stoichiometries using the endmember chemical formulae. reactions_from_formulae() is defined in burnman.tools.chemistry. Because there is only one possible equation I could pass reaction_stoichiometry[0] to equilibrium_temperature().
The code I suggested above should give you a temperature of 1556 K.

Best wishes,
Bob

10

Dear Gabriele,

I’ve now implemented the new functionality in burnman. You’ll have to install the bleeding edge version on your machine: e.g.

pip install git+git@github.com:geodynamics/burnman.git

or using git and a local version of the master branch.

Here’s a complete working code that solves both of your problems:

import burnman
from burnman import minerals, equilibrate
from burnman.tools.chemistry import equilibrium_temperature
from burnman.tools.chemistry import reactions_from_formulae
from burnman.tools.polytope import simplify_composite_with_composition


# Calculate the equilibrium temperature of an endmember reaction
fo = minerals.HP_2011_ds62.fo()
di = minerals.HP_2011_ds62.di()
en = minerals.HP_2011_ds62.en()
dol = minerals.HP_2011_ds62.dol()
CO2 = minerals.HP_2011_fluids.CO2()

assemblage = burnman.Composite([fo, di, en, dol, CO2])
reaction_stoichiometry = reactions_from_formulae(
    assemblage.endmember_formulae, assemblage.endmember_names, return_strings=False
)

P = 2.6e9
T = equilibrium_temperature(assemblage.phases, reaction_stoichiometry[0], P, 1500.0)
print(f"Equilibrium temperature of endmember reaction at {P/1.e9:.2f} GPa: {T:.0f} K")

# Calculate the equilibrium temperature CO2-in isopleth from assemblage
# endmembers: *en*, fs, fm, *odi*, *mgts*, cren, mess
opx = minerals.JH_2015.orthopyroxene([0.95, 0.0, 0.0, 0.04, 0.01, 0.0, 0.0])

# endmembers: *di*, cfs, cats, crdi, cess, jd, cen, cfm
cpx = minerals.JH_2015.clinopyroxene([0.95, 0.0, 0.01, 0.0, 0.0, 0.0, 0.04, 0.0])

# endmembers: *py*, alm, *gr*, andr, knor
gt = minerals.JH_2015.garnet([0.95, 0.0, 0.05, 0.0, 0.0])

# Create assemblage from the full solutions
assemblage = burnman.Composite(
    [fo, cpx, opx, dol, gt, CO2], [0.01, 0.01, 0.48, 0.48, 0.01, 0.01]
)

# Simplify the solutions because we are interested only in the CMAS system
assemblage = simplify_composite_with_composition(assemblage, assemblage.formula)

# Implement the equality constraints
equality_constraints = [["P", P], ["phase_fraction", [CO2, 0.0]]]

# Calculate the equilibrium from the equality constraints
sol = equilibrate(assemblage.formula, assemblage, equality_constraints)

print(f"Equilibrium temperature of CO2-in isopleth at {P/1.e9:.2f} GPa: {assemblage.temperature:.0f} K")

Best wishes,
Bob

11

Dear Bob,
what a professional work and reply!! Thank you very much.
I’ll work on it and we keep in touch.

Cheers,
Gabriele

12

Dear Bob first of all i hope everithing is going well and your 2024 started in the best way!

I came back after my honeymoon with my wife and for this reason i let in stand by my work with the equation on BurnMan.

I tried the script you furnished me and it works well for the first part:

import burnman
from burnman import minerals, equilibrate
from burnman.tools.chemistry import equilibrium_temperature
from burnman.tools.chemistry import reactions_from_formulae
from burnman.tools.polytope import simplify_composite_with_composition
fo = minerals.HP_2011_ds62.fo()
di = minerals.HP_2011_ds62.di()
en = minerals.HP_2011_ds62.en()
dol = minerals.HP_2011_ds62.dol()
CO2 = minerals.HP_2011_fluids.CO2()

assemblage = burnman.Composite([fo, di, en, dol, CO2])
reaction_stoichiometry = reactions_from_formulae(assemblage.endmember_formulae, assemblage.endmember_names, return_strings=False)

P = 2.6e9
T = equilibrium_temperature(assemblage.phases, reaction_stoichiometry[0], P, 1500.0)
print(f"Equilibrium temperature of endmember reaction at {P/1.e9:.2f} GPa: {T:.0f} K")

Then for the second part:

Calculate the equilibrium temperature CO2-in isopleth from assemblage

endmembers: en, fs, fm, odi, mgts, cren, mess

opx = minerals.JH_2015.orthopyroxene([0.95, 0.0, 0.0, 0.04, 0.01, 0.0, 0.0])

endmembers: di, cfs, cats, crdi, cess, jd, cen, cfm

cpx = minerals.JH_2015.clinopyroxene([0.95, 0.0, 0.01, 0.0, 0.0, 0.0, 0.04, 0.0])

endmembers: py, alm, gr, andr, knor

gt = minerals.JH_2015.garnet([0.95, 0.0, 0.05, 0.0, 0.0])

Create assemblage from the full solutions

assemblage = burnman.Composite([fo, cpx, opx, dol, gt, CO2], [0.01, 0.01, 0.48, 0.48, 0.01, 0.01])

Simplify the solutions because we are interested only in the CMAS system

assemblage = simplify_composite_with_composition(assemblage, assemblage.formula)

Implement the equality constraints

equality_constraints = [[“P”, P], [“phase_fraction”, [CO2, 0.0]]]

It is perfect until here.

Calculate the equilibrium from the equality constraints

sol = equilibrate(assemblage.formula, assemblage, equality_constraints)

After the script i found different errors here described:
AttributeError: ‘Solution’ object has no attribute ‘molar_fractions’

Is there any problem with the script?
Thank you in advance

13

I came back after my honeymoon with my wife

Congratulations on your marriage! I’m glad to hear you went on honeymoon with your wife, it would have been weird to go with someone else…

Is there any problem with the script?

Well, I can run it without problem…
I suspect you aren’t running it with the latest version of burnman. See the start of Message 10, above.

Best wishes,
Bob

14

Thanks Bob i solved the problem.
We keep in touch.
I think now we are able to prepare our work
Regards,
Gabriele

15

Dear Bob just some questions to complete our model.

if I wanted to calculate the molar proportions of each phase starting only from temperature and pressure, what script should I use?

I mean instead of creating the assemblage with the molar fractions (e.g. assemblage = burnman.Composite([fo, cpx, opx, dol, gt, CO2], [0.01, 0.01, 0.48, 0.48, 0.01, 0.01])), ask to burnman to report the molar fractions by inserting only P and T as equality constraints.

Everithing in order to create a dataset with the P-T couple and related expected molar fractions at that P-T equilibrium.

Thank you very much in advance.

16

Dear Gabriele,

That’s easy enough, just change:
equality_constraints = [["P", P], ["phase_fraction", [CO2, 0.0]]]

to
equality_constraints = [["P", P], ["T", T]]

And run equilibrate.

Then you can query the assemblage for the phase fractions and compositions.

Best wishes,
Bob

17

Dear Bob,
that’s what i tried, but something goes wrong.
This is exactly our run:

import burnman
from burnman import minerals, equilibrate
from burnman.tools.chemistry import equilibrium_temperature
from burnman.tools.chemistry import reactions_from_formulae
from burnman.tools.polytope import simplify_composite_with_composition

fo = minerals.HP_2011_ds62.fo()
di = minerals.HP_2011_ds62.di()
en = minerals.HP_2011_ds62.en()
dol = minerals.HP_2011_ds62.dol()
spl = minerals.HP_2011_ds62.sp()
mag = minerals.HP_2011_ds62.mag()
CO2 = minerals.HP_2011_fluids.CO2()

assemblage = burnman.Composite([fo, di, en, dol, CO2])
reaction_stoichiometry = reactions_from_formulae(assemblage.endmember_formulae, assemblage.endmember_names, return_strings=False)

P = 2.6e9
T = equilibrium_temperature(assemblage.phases, reaction_stoichiometry[0], P, 1500.0)
print(f"Equilibrium temperature of endmember reaction at {P/1.e9:.2f} GPa: {T:.0f} K")

opx = minerals.JH_2015.orthopyroxene([0.95, 0.0, 0.0, 0.04, 0.01, 0.0, 0.0])
cpx = minerals.JH_2015.clinopyroxene([0.95, 0.0, 0.01, 0.0, 0.0, 0.0, 0.04, 0.0])
gt = minerals.JH_2015.garnet([0.95, 0.0, 0.05, 0.0, 0.0])
assemblage = burnman.Composite([fo, cpx, opx, dol, gt, CO2], [0.01, 0.01, 0.48, 0.48, 0.01, 0.01])
assemblage = simplify_composite_with_composition(assemblage, assemblage.formula)
equality_constraints = [[“P”, P], [“T”, T]]
sol = equilibrate(assemblage.formula, assemblage, equality_constraints)

at this point when i try to run equilibrate, appears as follows:
LinAlgWarning: Diagonal number 13 is exactly zero. Singular matrix.
AssertionError: The highest upper bound for lambda is 1. (a full Newton step)

18

i also tried this way:
import burnman
from burnman import equilibrate
from burnman.minerals import SLB_2011

pressure = 1.8e9
temperature = 1400

composition = {‘Na’: 0.02, ‘Fe’: 0.2, ‘Mg’: 2.0, ‘Si’: 1.9,‘Ca’: 0.2, ‘Al’: 0.4, ‘O’: 6.81}
gt = SLB_2011.garnet()
ol = SLB_2011.mg_fe_olivine()
opx = SLB_2011.orthopyroxene()
cpx = SLB_2011.clinopyroxene()
assemblage = burnman.Composite(phases=[ol, opx, cpx, gt])
ol.set_composition([0.93, 0.07])
opx.set_composition([0.8, 0.1, 0.05, 0.05])
gt.set_composition([0.8, 0.1, 0.05, 0.03, 0.02])
cpx.set_composition([0.8, 0.05, 0.05, 0.05, 0.05])
equality_constraints = [(‘P’, 1.8e9), (‘T’, 1400.)]
sol, prm = equilibrate(composition, assemblage, equality_constraints)
print(f’It is {sol.success} that equilibrate was successful’)
print(sol.assemblage)

all right till now

But when i add dol and CO2 in the assemblage
burnman wants the term ‘C’ in composition
and then when i run equilibrate the response is:
“The bulk composition is not within the compositional space of the assemblage”.

I need to find the bulk composition for the complete assemblage (phases=[ol, opx, cpx, gt, dol, CO2])
{‘Na’: 0.02, ‘Fe’: 0.2, ‘Mg’: 2.0, ‘Si’: 1.9,‘Ca’: 0.2, ‘Al’: 0.4, ‘O’: 6.81}

19

Dear Gabriele,

Thanks for your detailed messages, they are very useful for understanding what you want to do.

This message:
LinAlgWarning: Diagonal number 13 is exactly zero. Singular matrix.
indicates that you are trying to solve a problem that has no solution.

For the system that you have set up, with endmembers for fo and dol, the line that we plotted last time is univariant. dolomite is stable on one side, and CO2 is stable on the other.

We can test this by appending the following three lines to the previous script:

equality_constraints = [["P", P], ["phase_fraction", [dol, 0.0]]]
sol = equilibrate(assemblage.formula, assemblage, equality_constraints)
print(f"Equilibrium temperature of dol-in isopleth at {P/1.e9:.2f} GPa: {assemblage.temperature:.0f} K")

You’ll see that the temperature is identical to the previous problem - dolomite becomes stable as CO2 becomes unstable.

So asking burnman to equilibrate this assemblage at any P-T point other than along the univariant line is meaningless.

Best wishes,
Bob