Magmasat degassing paths stop prematurely if one dissolved volatile goes to 0

[kaylai]

The reported issue

A user recently had an issue where the calculation of fully open-system degassing of a basanite resulted in the degassing path stopping when there was still ~1 wt% H2O dissolved in the melt. It seems that this occurred because the dissolved CO2 value became very small.

The final degassing step is at 170 bars, with dissolved H2O=0.985099 and dissolved CO2=2.474000e-17. Running a calculation of the saturation pressure with H2O=0.985099 and CO2=0 returns a consistent pressure at 170 bars.

Why this is happening

At the end of the code in the magmasat.py file, there are two lines that remove any rows from an open-system degassing data frame if CO2 or H2O in the liquid <=0. It looks like we added code to do this toward the beginning of VESIcal (the threshold at which to remove rows was updated from 0.000001, input at time of revisions of the manuscript, to 0 sometime shortly thereafter).

I don't remember why we made this decision. Presumably this was to stop crashing or otherwise erroneous behavior that we witnessed. However, it seems it is reducing the functionality of the model, so this needs to be reconsidered.

Tagging @simonwmatthews here to see if he has better recall than I do on this issue.

Complete working code example

mysample = v.Sample({'SiO2':  43.0832,
                    'TiO2':   4.2067,
                    'Al2O3': 15.1986,
                    'Fe2O3':  1.8962,
                    'Cr2O3': 0.0215,
                    'FeO':   10.1509,
                    'MnO': 0.1862,
                    'MgO':   8.4707,
                    'CaO':   9.7237,
                    'Na2O':   4.4341,
                    'K2O':    1.7525,
                    'P2O5':   0.8236,
                    'CO2': 0.61329,
                    'H2O': 2.116843
                    })

closed_degas=v.calculate_degassing_path(sample=mysample,temperature=1180, pressure=5880).result
open_degas=v.calculate_degassing_path(sample=mysample,temperature=1180,pressure=5880,fractionate_vapor=1.0).result

Result

Final line of the closed_degas dataframe is as expected:	Pressure_bars	H2O_liq	CO2_liq	XH2O_fl	XCO2_fl	FluidProportion_wt
1.0	0.062224	0.000016	0.891203	0.108797	2.669656

The open_degas dataframe, however, stops at P=170 bars with almost 1 wt% H2O dissolved in the melt:	Pressure_bars	H2O_liq	CO2_liq	XH2O_fl	XCO2_fl	FluidProportion_wt
170.2	0.985099	2.474000e-17	1.000000	6.724350e-15	0.221063

[simonwmatthews]

I will have a look at this. I have made some changes to the methods in the dev version I have been working on, so I will see if the issue still exists there.

[simonwmatthews]

I have had a look into this, and thought I would post an update before I lose track of what I have done. The issue doesn't seem to be the final lines which remove rows with < 0 H2O or CO2 (commenting out those lines changes nothing). Instead it seems to be that thermoengine is reporting 0 wt% fluid on these steps and so it never enters the if fl_mass > 0: loop.

I am now working on figuring out why it is doing that... will report back soon...

[simonwmatthews]

I think I have found a solution. For some reason, having a very small non-zero mass of CO2 in the bulk composition is causing thermoengine to not exsolve a fluid phase.

Adding this line in following the CO2 concentration adjustment seems to fix the problem for this case:

if _sample_dict["CO2"] < 1e-10:
    _sample_dict["CO2"] = 0.0

I will run the test suite etc then submit a pull request for this fix.

[kaylai]

Great news, thanks you!