[planetmariana]

Description

Adds planet.fO2_source = "from_mantle_redox". Instead of a fixed offset, the atmospheric fO₂ now comes from the melt's Fe³⁺/Fe²⁺ ratio, which is tracked as the mantle crystallizes. This follows Schaefer et al. (2024, JGR Planets) and uses the Deng et al. (2020) melt equation of state.

Potentially closes #653

New module src/proteus/interior_chem/

redox.py: update_melt_redox, called once per coupling step from interior_energetics/wrapper.py after the interior solve.

• Crystallization. Partitions Fe³⁺ and Fe²⁺ into the newly formed solid, using fixed partition coefficients (bridgmanite at P ≥ 22 GPa, Cpx/Opx mean below). The melt stays well mixed: one Fe³⁺/Fe²⁺ ratio for the whole melt.
• Metal saturation. After each crystallization step it checks 3FeO = 2FeO₁.₅ + Fe(metal) and brings any supersaturated melt back to a_Fe = 1. Metal is never redissolved.
• Surface fO₂. Hirschmann (2022) Eq. 21 (Schaefer Eq. 13) at 1 bar with the pressure term set to zero, as in Schaefer's fO2lowP_H22.m. It's evaluated at T_out = max(T_magma, outgas.T_floor), the temperature CALLIOPE and atmodeller solve at.
• ΔIW. Relative to tHirschmann (2021) IW buffer, also at T_out. Written to fO2_shift_IW_mantle, which CALLIOPE and atmodeller use as their fO₂ constraint.
• Solid threshold. Cells with melt fraction φ < 0.05 count as solid in every step (PHI_SOLID). The solver's own φ is not modified.
• Diagnostic fO₂ profile. For every melt cell, Eq. 13 at that cell's T and P, including the oxidation-reaction pressure term. It is stored but does not feed the outgassing.

disproportionation.py: equilibrium constant, metal activity and the solve for how much metal forms.

eos_deng.py: tabulated pressure integrals from the Deng et al. (2020) EOS, for both the disproportionation and the FeO → FeO₁.₅ oxidation reactions. Valid for T ≤ 5000 K (Deng's stated modelling range) and P ≤ 136 GPa.

Outputs

Helpfile columns: fO2_shift_IW_mantle, ferric_frac_mantle, a_fe_max_mantle, a_fe_max_cell_mantle, n_fe_metal_mantle, n_fe_metal_step_mantle, fe_metal_kg_mantle. The six besides fO2_shift_IW_mantle are diagnostics; fO2_shift_IW_mantle can be filled with zeros when resuming an older helpfile.

Interior snapshots:

• Aragog: log10_fO2_s and fO2_top_index are appended to data/<time>_int.nc, and also kept when the final snapshot is rewritten.
• SPIDER: a standalone data/<time>_redox.nc.

Config

• planet.ferric_fraction_initial (default 0.1).
• from_mantle_redox is rejected unless interior_energetics.module is spider or aragog.

Outgassing: outgas/wrapper.py, calliope.py and atmodeller.py now accept the new fO₂ source.

Docs: regenerated docs/Reference/config/ and docs/Reference/output_schema.json, and updated input/all_options.toml.

Validation of changes

Unit tests: tests/interior_chem/ (67 tests: redox, disproportionation, Deng EOS) and the new time-step cap tests in tests/interior_energetics/test_timestep.py. They cover:

• conservation of iron and oxygen in the metal step;
• that a supersaturated melt ends exactly at a_Fe = 1;
• an analytic low-pressure limit of the Eq. 13 pressure term (constant ΔV);
• the 1-bar / T_out rule for ΔIW, including a regression test that freezing the shallow melt doesn't shift ΔIW;
• the φ < 0.05 threshold: no melt iron in such cells, they're skipped by the metal check, and a cell crossing it crystallizes its remaining melt;
• that the results don't depend on whether cells are ordered surface-first (SPIDER) or core-first (Aragog);
• writing and reading the profile in .nc files;
• the a_Fe step cap: scaling to the target, inactive when a_Fe = 0, held at the minimum step, disabled at 0, decreases capped like increases.

Local test run (macOS 13.4, Python 3.12.4):

• ✅ Passed: tests/interior_chem, tests/utils/test_coupler.py, tests/interior_energetics/test_timestep*.py, tests/config, tests/tools/test_generate_config_reference.py. 607 passed.
• ⏭️ Not run: tests/interior_energetics/test_aragog.py, which fails 31 tests locally because of an outdated fwl-io. Clean main fails the same 31 in this environment.
• tools/generate_config_reference.py --check passes.

Cross-checks:

• The checked worked values come from the code's own functions: ΔG°, K and a_Fe for the metal reaction, the critical ferric fraction (f_crit = 0.069 at 3000 K and 100 GPa), and the surface fO₂.
• The tabulated pressure integral was compared with direct integration: median error 0.3–0.9%, at most 8.5% near 2 GPa.

Simulation runs:

Note

Merge fix included: while merging main, a return out on Aragog's success path was lost, and the redox helpfile keys were duplicated. Every accepted step was then counted as a failure, and runs stopped at the first step (T_surf = 0.0 K). Fixed in 696963a; the resulting aragog.py differs from main only by the snapshot hook.

Known limitations

• The metal check and the diagnostic profile skip melt hotter than 5000 K.
• The redox state isn't saved for resume: a resumed run restarts from ferric_fraction_initial.
• atmodeller applies the offset to its own IW buffer, which differs from O'Neill & Eggins by up to about 1 log unit at magma-ocean temperatures.

Checklist

• I have followed the contributing guidelines
• My code follows the style guidelines of this project
• I have performed a self-review of my code
• My changes generate no new warnings or errors
• I have checked that the tests still pass on my computer (with the environment caveats above)
• I have updated the docs, as appropriate (config and output references regenerated)
• I have added tests for these changes, as appropriate
• I have checked that all dependencies have been updated, as required (no new dependencies; netCDF4 and scipy are already required)

[codecov]

Codecov Report

❌ Patch coverage is 99.23896% with 5 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.75%. Comparing base (8849aea) to head (8194f19).

Files with missing lines	Patch %	Lines
src/proteus/interior_chem/redox.py	99.25%	0 Missing and 3 partials ⚠️
src/proteus/interior_chem/disproportionation.py	97.01%	1 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #948      +/-   ##
==========================================
+ Coverage   94.44%   94.75%   +0.31%     
==========================================
  Files         139      143       +4     
  Lines       23569    24217     +648     
  Branches     3944     4021      +77     
==========================================
+ Hits        22259    22947     +688     
+ Misses        882      879       -3     
+ Partials      428      391      -37     

Flag	Coverage Δ
unit-tests	93.53% <99.23%> (+0.15%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[timlichtenberg]

Thanks Mariana, this is a lot of careful work, and the supplement makes the model easy to follow. I ran the branch on the tutorial Earth case with a fixed-fO2 control and compared it with Schaefer et al. (2024). Before it leaves draft:

Use the Fischer et al. (2011) IW buffer, our default. The offset is taken against the O'Neill and Eggins buffer, but CALLIOPE adds it to its default Fischer et al. (2011) buffer, so the atmosphere gets a different fO2 than the tracker computes: +1.6 log units at 4240 K, +1.1 at 3000 K. Please compute the offset against Fischer et al. (2011) and test the value CALLIOPE receives.

Leave metal saturation out of this PR. With it on, Fe3+/FeT goes from 0.10 to 0.53 and the offset to IW+7, the whole-mantle outcome Schaefer et al. call implausible (p. 14), and it follows from the metal-forming cell sitting at the 5000 K validity limit. Let's focus this PR on the redox evolution without metal saturation, and come back to it in a follow-up PR once the EoS is extended far enough to allow metal saturation at the base.

A resume restarts the tracker. After a resume at 1722 yr, Fe3+/FeT goes from 0.155 to 0.03 for one step. The state needs to go into the snapshot.

Eq. 21 is also used far above its 1500 to 1900 K calibration, and the partition coefficients differ from the Schaefer et al. values above 15 GPa; a stated validity range and a benchmark against their mean curve would cover both. Please also remove proteus_update_20260929_164059_24312.log (it holds your home paths) and merge main. I can send the configs and the comparison plots.

[timlichtenberg]

Thanks Mariana, the resume works now. The crystallisation bookkeeping also holds up against Schaefer et al. (2024) when I drive update_melt_redox with their layer sequence and coefficients.

On the comparison, keep the melt-mass axis. PROTEUS has no magma-ocean base that moves up: the whole mantle is a mush that freezes at all depths together. The melt mass fraction that remains is the one axis both models share. Label it as that, with their base pressure as a second axis.

The disagreement comes from D(Fe2+). 0.85 is their bulk FeO coefficient for the fertile upper mantle (p. 5). In the lower mantle and the transition zone their Fe/Mg exchange gives a smaller D(Fe2+) than D(Fe3+) (p. 11), about 0.46 against 0.84 for the bridgmanite layers in their archived output, so their Fe3+/FeT falls there while ours rises. With their bulk coefficients per layer, the coupled run stays within 0.012 of their mean down to 10 % melt. I would split the validation into a benchmark in their set-up, as a unit test, and the coupled run with the offset compared at their adiabat temperature.

The coupled run with your coefficients (orange) and with theirs (green) against their 1000 draws; (c) and (d) show the two temperature conventions and how the PROTEUS mantle freezes.

Make both sides use Hirschmann (2021). The offset is against it now, but CALLIOPE still adds the offset to Fischer et al. (2011): the atmosphere is +0.35 log units at 4263 K to -0.65 at 700 K off the tracker's fO2. CALLIOPE#39 will make Hirschmann the default there (see my comment on it); this PR then needs that CALLIOPE version as its floor, and a test of the value CALLIOPE receives.

Metal saturation stays off by default. As we discussed, the code can stay in. We work on it in a follow-up PR. With it on, my run reaches 98 bar of O2 in the first 12 yr and ends at Fe3+/FeT 0.79. Please also rename the PR to match, without the metal saturation in the title.

Smaller: a cell that remelts returns no iron to the melt (Fe3+/FeT 0.219 against 0.192 with the return in my run), the page should state the 1468 to 2197 K calibration of Eq. 21, and redox.py.cryst can go. The PR description and the redox.py docstring still describe the earlier state (the O'Neill and Eggins buffer, the metal step always on, tools/redox_step.py), and 26 comment blocks are over 3 lines. I can send the configs and the benchmark script.

[planetmariana]

Thanks for the review @timlichtenberg these are the latest changes from your comments. I'm pinning it here now that I'm done, so please take your time and have a look once you're back in the office:

• Hirschmann on both sides: Is done from this PR on PROTEUS side, but needs first merge CALLIOPE#39
• Metal saturation: off by default, code kept for the follow-up.
• Remelting: a remelting cell now returns iron from the solid it formed, and a remelt after full solidification restarts the tracker. Tests cover conservation and the round trip to f0.
• Smaller points: Eq. 21 calibration range added to the page, redox.py.cryst removed, docstring updated, no added comment block over 3 lines.

Now regarding why PROTEUS was bending down the Fe3+/FeT ratio, based on your previous comment, I found out this:

Difference from Schaefer et al. (2024) late in crystallisation. PROTEUS keeps Fe3+/FeT bending down at the end, while S24 bends up. This comes from how solid forms, not from the partition coefficients:

• Tracker on S24's sequence: fed S24's own crystallisation sequence (their Mr_r_P_wholeEarth.mat, bottom-up), the tracker gives S24's dip and late rise with both coefficient sets.
• S24: the late solid forms only at shallow pressure, in olivine-rich layers with D(Fe3+) − D(Fe2+) < 0, and the trapped liquid of each layer is closed off. Fe3+ stays in the melt, so the ratio rises.
• PROTEUS: the melt fraction comes from the energetics, and all melt, including deep melt, stays as one reservoir. About half of the late solid still forms at depth with bridgmanite-like D, so the effective D(Fe3+) − D(Fe2+) stays positive and the ratio keeps falling.
• What decides it: treating cells below a critical melt fraction as solid brings back the S24-like rise, so it comes down to whether deep mush stays in contact with the ocean. S24 and PROTEUS take the two end-members.
• Coefficients: without trapped liquid, as in PROTEUS, mgsio3 follows the S24 mean more closely than bulk_BSE on S24's sequence (RMS 0.005 against 0.011 in Fe3+/FeT), so mgsio3 stays the default.

[timlichtenberg]

Thanks Mariana, the remelting now returns its iron and the resume keeps the tracker, and the default path is unchanged. Two things before this can go in, and a few smaller ones.

CALLIOPE still adds the offset to Fischer et al. (2011). The tracker computes the offset against Hirschmann (2021), so the atmosphere sits -0.65 log units off the tracker's fO2 at 700 K and +0.35 at 4263 K. The current test only checks that the offset arrives in CALLIOPE, not the fO2 it then uses.

How I would fix it: wait for CALLIOPE#39 to be released, then raise the fwl-calliope floor in pyproject.toml to that release in this PR. Add a test that runs the CALLIOPE call with a known tracker state and compares the absolute log10 fO2 CALLIOPE uses with the tracker's at two temperatures, one cold and one hot. Once CALLIOPE#39 is merged, make the buffer a PROTEUS setting that the tracker and CALLIOPE both read (fO2_model on the CALLIOPE side), so the two cannot drift apart, and refuse a config where they would differ.

Late in the freeze the tracked melt holds more iron oxide than melt mass. In a run to a melt fraction of 1e-4, FeO over melt mass goes from 0.09 at the start to 0.70 at the default phi_crit of 0.01 and above 1 below a melt fraction of 0.005. Over the same stage the offset passed to CALLIOPE falls from +3.1 to -1.4, so this state sets the outgassing.

The cause is the explicit step: with a constant D below 1 it has no bound. It can also drive n_fe3_melt negative with the bulk_BSE coefficients when a step freezes most of a layer, and it leaves iron in a melt of zero mass when the last melt freezes.

How I would fix it: replace the explicit update by the per-step closed form, either batch, n(1-s)/(1-s+sD), or Rayleigh, n(1-s)^D. Both stay within [0, n] for any step and empty the melt at s = 1, so the result no longer depends on how coarse the late steps are. Then stop updating the tracker below a stated melt fraction and keep the last state, and add tests for a step that freezes the whole layer and for the melt FeO staying below the melt mass.

Smaller:

• The PR description still describes the earlier state (metal check after every step, PHI_SOLID 0.05, resume from ferric_fraction_initial, O'Neill and Eggins buffer, 5000 K cut).
• _reject_reserved_fO2_source (_planet.py:184) is never called, and its message says the source has no runtime path yet; it can go.
• The redox pages carry dated history (docs/Explanations/redox.md:223, :454 and docs/Validation/interior_chem/redox.md:7, :48) and a figure that has not been regenerated; redox.md:49 says metal_saturation must be false, while the config accepts true.
• Two tests cannot fail: test_reaction_gibbs_is_insensitive_to_the_reference_state_choice builds shifted from direct plus and minus the same term, and test_planet_validators.py:76 asserts on its own parametrize list.