Total surface pressure being reset to 1 bar when allowing partial pressures to update from helpfile

[nichollsh]

Currently in the master branch, the flags passed to SPIDER which set the initial partial pressures of each volatile *_initial_atmos_pressure are not updated at runtime. They are left at the initial values provided in the configuration file, and as a result the partial pressures don't really change throughout the simulation. They do change slowly, as they are allowed to evolve during each SPIDER call, but then are reset at the start of the next SPIDER call. This is evident in the plot_global graph (see figure below, specifically subplot D).

If we introduce a patch to the RunSPIDER() function which updates the value passed via this flag according to the helpfile, we then see volatile partial pressure evolution over time. This is included in a branch called patch_230512_pressure. The results of running the code with this patch are shown in the figure below, with the same initial conditions as the one above. Allowing the partial pressures to evolve in this way also affects the solidification time, which I think physically makes sense.

However, this patch introduces a new issue where the total pressure is being reset to 1 bar after the init steps have been completed, with the partial pressures being scaled such that mixing ratios are conserved. I don't know what's causing the total pressure to be reset to 1 bar, but it's something occurring either in SPIDER itself, or in the way its JSON files are being read.

I've also included an Excel spreadsheet showing the evolution of both COUPLER_options and the helpfile, for both of these simulations: evolution_pathways.xlsx

This issue exists on both MacOS and Linux installations of the code.

[timlichtenberg]

My first guess is that changing _initial_atmos_pressure during runtime violates volatile mass balance within SPIDER. @djbower: Can you change _initial_atmos_pressure in subsequent SPIDER calls, when restarting SPIDER from a new data dump? Let's say H2O_initial_atmos_pressure = 1.00e5 at the t = 0 (fully molten); does this parameter have to be conserved when restarting from a later data dump and no volatile loss mechanisms (such as escape) are active?

[djbower]

I will check next week how the coupling is actually performed in practice. Formally of course, pressure is never conserved, but rather elemental mass/moles is conserved (e.g. total H, which can be partitioned between H2, H2O, CH4 etc.). Either way, SPIDER has multiple ways of restarting the volatiles, so I expect a simple adjustment somewhere in the workflow will fix the problem.

On Fri, 12 May 2023, 15:29 Tim Lichtenberg, @.***> wrote:

My first guess is that changing _initial_atmos_pressure during runtime violates volatile mass balance within SPIDER. @djbower https://github.com/djbower: Can you change _initial_atmos_pressure in subsequent SPIDER calls, when restarting SPIDER from a new data dump? Let's say H2O_initial_atmos_pressure = 1.00e5 at the t = 0 (fully molten); does this parameter have to be conserved when restarting from a later data dump and no volatile loss mechanisms (such as escape) are active?

— Reply to this email directly, view it on GitHub https://github.com/FormingWorlds/PROTEUS/issues/37#issuecomment-1545748197, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOYQUYJ6NJSHALMD3QNOVDXFY3KXANCNFSM6AAAAAAX7G6QW4 . You are receiving this because you were mentioned.Message ID: @.***>

[djbower]

In SPIDER, *_initial_atmos_pressure is only used during the construction of the initial condition for the atmosphere if IC_ATMOSPHERE=3. Then the code simply uses the predefined partial pressures to back-compute all other relevant quantities. If IC_ATMOSPHERE is not 3, then *_initial_atmos_pressure is conformed to be consistent with however the initial condition was specified. Once the partial pressures of all species have been determined, it is this quantity that is evolved during the time-stepping (*_initial_atmos_pressure remains constant because it is not used).

Hence, if you are coupling the species to SPIDER using the partial pressure, indeed you must update *_initial_atmos_pressure and set IC_ATMOSPHERE=3 every time you pass control of the simulation back to SPIDER.

"However, this patch introduces a new issue where the total pressure is being reset to 1 bar after the init steps have been completed"

I can't immediately see this as a SPIDER issue. Check the individual JSON files (particularly the first output) to confirm that the partial pressures you are putting in are the same as those that appear in the JSON. The total atmosphere pressure, which should be equal to the sum of the partial pressures, is also given in the field atm_struct_pressure. Also, the volatile masses are output, so you can take a look at those. In particular, see if total volatile mass is conserved in your combined simulation where you are regularly restarting SPIDER. And if it isn't conserved, let's see how much mass is lost or gained for each species considered.

[timlichtenberg]

Thanks! That makes sense, but @nichollsh and I discussed about it, and we are not sure how SPIDER is handling *_initial_atmos_pressure with IC_ATMOSPHERE=3 in the case of a restart at, say, t = 10; not at the very start. Does it find the matching volatile mass in the interior/melt with *_initial_atmos_pressure be treated as fixed, or does it treat *_initial_atmos_pressure as total volatile mass, and redistributes it dynamically? When doing that, and restarting PROTEUS with *_initial_atmos_pressure informed by the last time from the previous SPIDER call, we seem to lose volatile mass, see the plot by Harrison below.

[nichollsh]

[djbower]

Whenever SPIDER is started with IC_ATMOSPHERE=3 the value of _initial_atmos_pressure is used to determine the total mass of the volatile in the system. This same operation happens regardless of the restart time (there is no back-determination to calculate what the initial would have been at time zero).

When the volatile mass is computed from _initial_atmos_pressure, it effectively computes the volatile mass in each of the three reservoirs: the interior melt, the interior solid, and the atmosphere. Then it sums all the reservoirs to give the total. This total is then conserved for the subsequent evolution until SPIDER is halted. But once it is restarted, it will again recompute the total volatile mass if IC_ATMOSPHERE=3.

[timlichtenberg]

Ok, if understand this correctly (e.g. for H):

H_tot = H_melt + H_solid + H_atm,

where H_atm is set by _initial_atmos_pressure. This means, the code does (in principle)

1. H_atm = _initial_atmos_pressure
2. H_melt = H_atm * {Atmosphere-melt volatile partition coefficient}
3. H_solid = H_melt * {Solid-melt volatile partition coefficient}
4. H_tot = H_melt + H_solid + H_atm

So, if _initial_atmos_pressure is set to the value of the partial pressures outputted by SPIDER from the last data file, then it should end up with the same total overall mass. Is that overall logic correct?

[djbower]

Yes that's correct. Not only will _initial_atmos_pressure need to be the same, but also the global melt fraction since this is used to determine the mass of melt and solid. Of course, please run a quick check and confirm this is the case.

[timlichtenberg]

Ok, thanks, need to do some more tests.