SPIDER call redox chemistry

[timlichtenberg]

Goal: Implement SPIDER call sequence to specific redox state + H/C ratio instead of total volatile budget or partial pressure

@djbower Could you provide us with an example to call sequences to SPIDER that do that so we can implement this into PROTEUS?

Are there incompatibilities we should be aware of with other call variables? The list of used variables is here (line 67 ff.): https://github.com/FormingWorlds/PROTEUS/blob/master/init_coupler.cfg

[djbower]

See parameters.c around line 1032.

  Ap->OXYGEN_FUGACITY = OXYGEN_FUGACITY_NONE;
  {
    PetscInt OXYGEN_FUGACITY = 0;
    PetscBool OXYGEN_FUGACITYset = PETSC_FALSE;
    ierr = PetscOptionsGetInt(NULL, NULL, "-OXYGEN_FUGACITY", &OXYGEN_FUGACITY, &OXYGEN_FUGACITYset);
    CHKERRQ(ierr);
    /* only set OXYGEN_FUGACITY if there are reactions, since fO2 is
       used to also adjust the total atmosphere pressure */
    if (OXYGEN_FUGACITYset && Ap->n_reactions)
      Ap->OXYGEN_FUGACITY = OXYGEN_FUGACITY;
  }

  /* default offset, assuming IW buffer is 0.5 from Sossi et al. (2020) */
  Ap->OXYGEN_FUGACITY_offset = 0.5;
  ierr = PetscOptionsGetScalar(NULL, NULL, "-OXYGEN_FUGACITY_offset", &Ap->OXYGEN_FUGACITY_offset, NULL);
  CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

Set OXYGEN_FUGACITY to an integer according to the models specified in atmosphere.c::set_oxygen_fugacity. You can also specify OXYGEN_FUGACITY_offset to give the log10 offset relative to the buffer. If you generally search in the code for OXYGEN_FUGACITY you will find the required functionality.

Currently there is no way in SPIDER to specify a H/C ratio directly. You have to determine either the pressures of the species or the total mass (kg) and then feed those into the code. In which case, the coupler will need to do the conversions.

[timlichtenberg]

Awesome, thank you, that is even better! So we can feed the H-C-N-O compounds via H2O_initial_total_abundance (these are ppm wt though?) and set OXYGEN_FUGACITY or OXYGEN_FUGACITY_offset. I am not sure I understand the difference between the latter two parameters though.

[djbower]

OXYGEN_FUGACITY sets the choice of the buffer (e.g. Iron-Wustite or whatever) whereas OXYGEN_FUGACITY_offset is the log10 offset relative to the chosen buffer. So for example you can choose the IW buffer and then set the offset to '-1' in order to set the fugacity at IW-1.

Yes initial total abundance is ppm by weight, which is the mass of the volatile divided by the total mantle mass (which could be a combination of melt or solid).

[timlichtenberg]

Great, thank you! :-)

[nichollsh]

Hi @djbower , this approach seems to work when IC_ATMOSPHERE=3 and volatile partial pressures are set. However, the model crashes when using IC_ATMOSPHERE=1 and the total abundances (ppmw) are set.

The error occurs when the model tries to solve for the partial pressures from the total abundances. I've provided more information below and attached a couple of options files that have produced this error (file extension changed to txt, as GitHub didn't allow opts). Could you tell me what I might be doing wrong?

Run command: ./spider -options_file test_ic_abundance.opts

Error message: Matrix type mffd does not have a multiply transpose defined or is symmetric and does not have a multiply defined.

Error sequence:

[0]PETSC ERROR: #1 MatMultTranspose() at /home/n/nichollsh/PROTEUS/petsc/src/mat/interface/matrix.c:2625
[0]PETSC ERROR: #2 SNESSolve_NEWTONTR() at /home/n/nichollsh/PROTEUS/petsc/src/snes/impls/tr/tr.c:429
[0]PETSC ERROR: #3 SNESSolve() at /home/n/nichollsh/PROTEUS/petsc/src/snes/interface/snes.c:4666
[0]PETSC ERROR: #4 solve_for_initial_partial_pressure() at /home/n/nichollsh/PROTEUS/SPIDER/ic.c:940
[0]PETSC ERROR: #5 set_ic_atmosphere_from_initial_total_abundance() at /home/n/nichollsh/PROTEUS/SPIDER/ic.c:828
[0]PETSC ERROR: #6 set_ic_atmosphere() at /home/n/nichollsh/PROTEUS/SPIDER/ic.c:468
[0]PETSC ERROR: #7 set_initial_condition() at /home/n/nichollsh/PROTEUS/SPIDER/ic.c:54
[0]PETSC ERROR: #8 main() at /home/n/nichollsh/PROTEUS/SPIDER/main.c:78

test_ic_abundance_simple.txt test_ic_abundance.txt

[djbower]

It is often difficult for the solver to find an initial partial pressure from total abundances if the total abundances are far from the basin of convergence. Unfortunately, there's not really a good way to deal with this, other than to trial-and-improvement the initial condition.

The approach that I previously took was to start with just H2O (in which case you don't need OXYGEN_FUGACITY or OXYGEN_FUGACITY_offset since you have no reactions and only one species). Then you make sure that the system can solve for the initial partial pressure. Once this works, you can add in H2. You know the ratio of pH2O/pH2 from equilibrium chemistry at a given temperature so you can try and enforce this ratio from the outset. For example, for a given gamma = pH2O/pH2 and pH2O, you can compute the total mass of H in the system using a simple python script. Then you can solve to get the initial H2 and H2O partial pressures.

Then you can add in the carbon system, since the carbon system will only change the partial pressures of the hydrogen species by a small amount. And then basically you keep iterating and adding more species to ensure you are close enough to the equilibrium state that the solver finds a solution. This limitation was ultimately why I switched to using the partial pressure for the initial condition, rather than trying to back-solve from the initial total abundance. Because practically, it was easier to evolve partial pressures from a known "correct" state than to ask SPIDER to solve for the initial partial pressures.

Another approach that is helpful is to consider what is the main background gas. In most cases at high temperature it is CO, in which case you can approximation mu_bar = mu_CO in C.6. of Bower et al. (2022). This allows you to decouple the different systems (e.g. C and H) and you can treat them as independent equation sets to get a rough solution using a simple python script.

So in short, yes there can be problems solving for the initial partial pressures of the volatiles, and this is legacy of the fact that SPIDER was not designed from the outset with a sophisticated chemical solver. But once you help SPIDER to be within the basin of convergence by presolving for the initial pressures (even roughly), SPIDER does a good job of integrating from that point throughout the model.

[djbower]

I've attached the python script I wrote to determine the initial partial pressures for SPIDER based on a C/H ratio, H mass, and fO2. This duplicates much of the code you see in SPIDER. Obviously you have to make sure that this script and Python are using the same solubility law, fO2, etc. There's also an option to include nitrogen. To the best of my knowledge this script should work out of the box. This will allow you to then feed the solved partial pressures into SPIDER directly and avoid asking SPIDER to back-solve.

spider_volatiles_solver.py.zip

[nichollsh]

That's fantastic - your explanation makes sense. Thanks a lot for your help, Dan!

Since I was having no problems with using partial pressures, I had already considered writing a short script to compute them from the above parameters. However, it seems that you've already done exactly that, which is very helpful. The script itself runs fine out of the box, and is reproducing the expected behaviours.

Would it be okay with you if this code was adapted into PROTEUS? It would be ideal to include these parameters inside the PROTEUS configuration file, instead of the partial pressures, as it makes for more meaningful comparisons.

[djbower]

Sure, feel free to include the code in PROTEUS or adapt and amend it as required. If you happen to be keeping a reference list for codes, perhaps you can add this paper:

https://www.sciencedirect.com/science/article/pii/S0012821X22005301

Simply because I wrote the script that I sent you for the above paper (notably the MCMC part).

[nichollsh]

Thanks - I'll make a note of the paper.