[IKisvardai]

Description

Adds solid-phase volatile trapping during mantle crystallisation. The trapped
melt fraction comes from the drainage of pore melt across the freezing front
that the Aragog interior resolves. Relates to #519.

outgas.trap_mode takes 'none' (default) or 'front'. With 'none' the
trapping step does not run, desiccation and the escape floor act as on main,
and the per-element closure check is not applied. One fix applies with
trapping off as well: a frozen mantle's atmosphere is now debited by what
escape took from each element total (see Bookkeeping), where main left
1.8e-5 of every step's loss in the atmosphere with CALLIOPE. The golden
trajectory differs from main only by new constant-zero columns, and a
helpfile written by main resumes on this branch.

Physics

• Mass balance. Each crystallisation step moves
  [(1 - F_tl) D_Z + F_tl] C_Z dM_RM of every dissolved species from the melt
  into the solid mantle (Sim, Hirschmann & Hier-Majumder 2024, Eq. 6). Here:

  • C_Z is the concentration from the previous chemistry solve, taken on the
    melt it dissolved into;
  • dM_RM = M_mantle [Phi(t-1) - Phi(t)];
  • D_Z comes from outgas.D_const_<species>. Water is 0.0017 and every
    other species is 0, so the others are buried through the interstitial melt
    alone. That includes the noble gases dissolved in the melt, at
    F_tl C_Z dM_RM.

  The step is capped at what the melt holds.

• Trapped melt fraction. The front runs from the node where the solver's
  melt fraction crosses interior_energetics.rfront_loc down to the porosity
  floor trap_phi_min. A parcel enters it at the top at porosity phi_top.
  Over the residence time t_res = L / v_f it follows
  dphi/dt = -phi / max(L / w_D(phi), tau_s), and it keeps phi(t_res),
  where:

  • w_D = (1 - phi) |drho| g F(phi) / eta_melt is the Darcy speed. F is
    Aragog's own three-regime permeability over porosity,
    aragog.eos.mobility_function (Bower et al. 2018), and the drag is the
    melt viscosity;
  • tau_s = mu_s / (|drho| g L) is the matrix deformation time, with the mush
    viscosity trap_mush_log10visc (default 1e20 Pa s);
  • v_f = |dM_RM| / (4 pi r_base^2 rho_s dt) is the front speed from the
    crystallised mass.

  t_res depends on the front speed, not on the step length. The retained
  volume fraction is converted to a mass fraction.

• Guards. A front the mesh cannot resolve takes the no-drainage bound,
  which is the entry porosity capped at rfront_loc, and logs the cause. The
  causes are:

  • a front too thin to resolve;
  • a mush split into separate layers;
  • melt denser than the solid inside the front;
  • a front reaching the surface.

  If the interior profiles are unavailable, the step buries at D_Z alone
  (F_tl = 0), on a branch of its own.

• Remelting. The solid that remelts first is the solid that froze last.
  A burial ledger (outgas/burial_ledger.py) records each step's burial
  against the interval of melt fraction the mantle froze through, and a step
  on which Phi_global rises to Phi returns what is recorded below Phi:

  • Freezing and remelting the same interval cancel exactly.
  • A mantle that remelts completely gets back everything trapped, however
    the melt changed in between.
  • The ledger is rebuilt from Phi_global and the _kg_trapped columns of
    the completed rows. It keeps no state of its own, survives a resume, and
    the rebuild takes about 35 ms per remelt step on a 5000-row history.
  • Trapped mass the history does not account for, as after resuming a
    helpfile written before trapping, returns in proportion to the solid that
    remelted, (Phi(t) - Phi(t-1)) / (1 - Phi(t-1)).

  Remelting steps report branch 7.

• Front profiles. The front is located and drained on Aragog's basic
  nodes, from the core-mantle boundary up, with the melt fraction, mixture
  density, phase-boundary densities and gravity that Aragog 26.10.03 writes
  there (phi_basic, rho_basic, rho_solid_b, rho_melt_b, g_b). No
  solver internals are read.

• Porosity. Inside the mush, porosity comes from the density lever rule
  on those densities with the sign of the contrast kept. That gives the exact
  volume fraction for Aragog's two-phase mixture, whichever phase is denser.
  Aragog's porosity_b is the same rule with a soft clip where the melt is
  lighter, but reads about 0 where the melt is denser, as its docstring warns,
  so it is not used for the drainage.
  Single-phase nodes take their melt fraction.

Bookkeeping

• Trapped mass is recorded in new {species}_kg_trapped and
  {element}_kg_trapped columns and counted in _kg_solid. There is one
  column for each volatile species, volatile element and noble gas. The rest
  of _kg_solid stays with the chemistry backend, for example atmodeller's
  graphite.
• Around each chemistry solve, the trapped mass is taken out of the element
  totals and the solid columns, then put back after the solve. A try/finally
  context manager does the swap, and the species totals are rebuilt from the
  new partition.
• A species is buried only when every one of its elements has a positive
  whole-planet total, so dust-level melt mass of an element the run does not
  carry stays where the chemistry put it.
• Escape and the desiccation gate see only the mass they can reach.
  O_kg_total includes the trapped oxygen under both fO2 sources.
• Desiccation follows main (Keep the escape floor off a frozen mantle and book desiccation #953): it empties the atmosphere, books the
  removed mass in M_desiccated, and keeps the solid share, trapped mass
  included. With trapping on it keeps the melt as well, so mass a remelt
  returns to a desiccated planet stays in the planet however many steps the
  remelt takes.
• After a remelt step on a desiccated planet, desiccation is checked again,
  and the flag clears once the mass held in the melt reaches mass_thresh.
• The escape floor reads the reachable mass, the total less what trapping
  locks in the solid, as the desiccation gate does, and empties it down to
  the locked mass. It is off on a frozen mantle, as on main, and on a step on
  which a remelt returned buried mass, so that a slow remelt keeps what a
  fast one keeps.
• On a frozen mantle, escape debits each element total by its share of the
  escaping mass, its atmospheric mass over the sum of the element masses in
  the atmosphere. run_crystallized now scales the atmosphere by the
  fraction of that same sum that remains, so each element's atmosphere loses
  exactly what escape took from its total. It used (M_atm - esc) / M_atm,
  and M_atm, the sum of the species masses, is 1.0000177 times the element
  sum with CALLIOPE's molar masses: 1.8e-5 of each step's loss stayed in the
  atmosphere, and as escape drew a total down towards what the solid holds
  that gap crossed the closure limit. Species columns and pressures scale by
  the same fraction, so the composition and the ratio of M_atm to the
  element sum are kept.
• With trapping on, assert_mass_conservation(check_element_closure=True)
  checks total = atm + liquid + solid for every element at every step. A
  mismatch is refused once it exceeds both 1 kg and the chemistry solver's
  relative tolerance times the total.
• A helpfile without the trapping columns resumes. The 37 trap_* and
  *_kg_trapped columns are read as zero through the same path as main's
  RESUMABLE_ZERO_FILL_KEYS: zero is their value in a run that never
  trapped. Any other missing column is still refused.

Other changes

• A config-load validator refuses trap_mode = "front" with any interior other
  than Aragog.
• New helpfile columns:
  • trap_* diagnostics: branch, fraction, timescales, front geometry and
    speed, and buried mass per step and cumulative. trap_front_courant
    reports the front's advance over one step, in units of its own thickness;
  • the _kg_trapped columns.
• The init stage never traps.
• Zalmoxis takes its dissolved-H2 target as total minus atmosphere minus
  solid.
• The output-reference tools understand the registry loops over vol_list,
  vol_element_list and noble_gases.
• The config and output references are regenerated, the options appear in
  input/all_options.toml, and there are validation pages for trapping.py
  and compaction.py.

Modelling choices, for review

1. w_D carries (1 - phi), which Aragog's relative_velocity omits. Both
   give the melt velocity relative to the matrix: McKenzie (1984, Eq. 3.7)
   takes the buoyancy against the mixture, (1 - phi) (rho_s - rho_f) g, and
   Abe (1995, Eq. 40) against the solid. The factor lengthens the percolation
   time by 1 / (1 - phi), 1.43 at phi = 0.3, and in the five test cases
   below raises the retained fraction by at most 0.0095. McKenzie's solution
   assumes touching grains, so melt fractions below about 0.2, while the front
   is entered at rfront_loc, 0.4 by default; Measure the (1 - phi) Darcy factor in the gravitational separation velocity aragog#67 tracks
   the same factor for Aragog.
2. There is no front-thickness limit by default
   (trap_max_front_fraction = 1). Where percolation controls, the retained
   fraction depends on w_D / v_f and not on L. A value below 1 sends
   thicker fronts to the no-drainage bound.
3. The front speed comes from the crystallised mass dM_RM, the same increment
   the flux uses, rather than from M_mantle_solid. The two therefore cannot
   disagree about whether a step froze.
4. The drag viscosity is always the melt viscosity, whatever the solver's
   separation_viscosity.
5. The scheme is global: one front per step and no per-node state. The
   burial ledger is kept in melt-fraction coordinates, not radius: it orders
   the solid by when it froze, which is what decides what remelts first in a
   global model.
6. The closure check's absolute tolerance is a fixed 1 kg. It admits the
   trace melt mass the chemistry leaves of an element the run does not carry
   (6.7e-2 kg in the Earth run with He and Ar) and is independent of
   mass_thresh, so a break below that setting is still caught.
7. With trapping on, a desiccated planet keeps its melt as well as its solid;
   the result of a remelt then does not depend on how many steps it takes.
8. The escape floor is off on a step on which a remelt returned buried mass:
   the element is being resupplied, and flooring the small release each step
   would delete it whenever it is below mass_thresh.
9. With freeze_volatiles, burial goes on between the freeze and
   Phi_global = 0. The last melt crystallises and buries at F_tl while the
   atmosphere stays out of the exchange. What the burial leaves in the melt
   reservoir stays there, as the whole reservoir does on main.

Known limitations

• The front scheme assumes a thin front. When the mush spans most of the
  mantle, v_f is only an effective speed, and with the dense-melt and
  surface-node guards now able to fire, such steps take the no-drainage bound.
• The PALEOS phase-boundary tables make the melt denser than the solid below
  about 5 GPa. That looks like a table artifact, but trapping takes it as
  given (the dense-melt guard).

Validation of changes

Unit tests (tests/outgas/test_trapping.py, test_compaction.py,
test_wrapper.py, tests/utils/test_coupler.py, tests/config,
tests/interior_struct/test_zalmoxis.py):

• Eq. 6 is pinned against hand-computed values and its analytic limits, and
  the supply cap is tested through run_trapping.

• Remelting, tested on a mantle frozen from Phi = 1 to 0.5, which buries
  4.333e19 kg of water, then degassed a hundredfold:

  • remelting to 0.75 returns half of every trapped species and element;
  • remelting to 1 returns all of it, leaving the trapped and solid reservoirs
    at zero;
  • sizing the release from the degassed melt would have returned 1.98% and
    left 98% trapped;
  • a mantle remelting from fully solid (Phi = 0, no melt) still releases
    its share.

• With that one-row history the ledger has nothing to replay, so these
  releases test the fallback share; helium runs through the same partial and
  complete remelt, coming back once, with its moles following the masses
  after the withhold and restore.

• A closed melt frozen from Phi = 0.6 to 0.4 and cycled a hundred times
  between 0.41 and 0.40 keeps its trapped water to 1e-12. Releasing the mean
  of the solid instead raises it by 4.4% a cycle and triples it in a hundred.
  The ledger itself is checked against an independent fine-grid ledger (2e5
  melt-fraction bins) on a seeded random history of 400 steps.

• Dissolved helium is buried at F_tl C_Z dM_RM, withheld from the chemistry
  and restored with its moles.

• With a zero N budget, 6.7e-2 kg of N2 and 1e-2 kg of NH3 in the melt are
  not buried while water still is. The closure check passes the dust-level N
  row and refuses it with no absolute tolerance, and refuses a 2 kg break on
  a dust total, a 1e15 kg break on a 5e14 kg total (below a mass_thresh of
  1e16 kg) and a 10% break on a planetary one.

• A desiccated planet with H, O and Ar trapped that remelts half its solid
  is no longer desiccated, with every element total kept, and kept also
  while it stays flagged.

• A desiccated planet with 4e18 kg of H and 3.2e19 kg of O trapped, remelting
  from Phi = 0.30 to 0.50 in 1000 steps of 2e-4 that each release less than
  mass_thresh, keeps every total and ends with 5/7 of its water trapped and
  1.03e19 kg in the melt, as 4 steps of 0.05 do.

• The same slow remelt runs through Proteus.start on the dummy
  configuration: 16 steps of 0.0125, each releasing less than mass_thresh,
  with the flag set over the first remelt step and the released hydrogen
  below the threshold for the rest of the run. The totals hold on every row
  and the run ends with the trapped 5/7 and the outgassed 2/7 that 4 steps of
  0.05 reach. The planet also holds 1.5e18 kg of trapped N2, whose released
  2/7 stays below mass_thresh: it is kept while the remelt lasts and, once
  it stops, floored to the locked 5/7, which a whole-total floor would not do.

• A frozen H2O-N2 atmosphere whose species masses are 1.0000177 times its
  element masses, as with CALLIOPE, goes through 60 alternating escape and
  frozen steps of a quarter of the atmosphere each: every element closes to
  1e-10 while the N total falls to what the solid holds. Through
  Proteus.start, the dummy mantle is frozen with freeze_volatiles and
  trapping on and escapes for 20 steps at that ratio, with every element
  closing to 1e-10 on every row; with the old scaling the loop stops on the
  closure check for H at 1.07e-4.

• Desiccation with trapping on keeps the solid and the melt and books the
  atmosphere in M_desiccated; with trapping off it empties the melt as on
  main.

• Two molten-mantle runs that differ only in the trapped N, 5e14 or 1.1e16
  kg, both empty a reachable 9e15 kg below mass_thresh to the locked mass,
  where a whole-total floor would empty one and keep the other. On a remelt
  step the floor is off; a step that released nothing keeps it.

• The trapping step's lever-rule porosity on Aragog's densities equals
  Aragog's own porosity_b formula to 2e-5 where the melt is lighter, and
  keeps the melt volume where Aragog's reads 0.

• The same front crossed at the same speed drains the same F_tl over a step
  five times longer, with the front advancing 1.99 of its own thickness.

• A helpfile written without the 37 trapping columns loads with them zeroed.
  A missing core column is still refused.

• Per-element closure holds after burial and after release.

• run_outgassing is driven with a stand-in backend. It hands the chemistry
  the reachable inventory, keeps atmodeller-style graphite separate from
  trapped mass, rebuilds species totals, and restores the trapped mass when
  the solve raises.

• The permeability is checked against Aragog's own
  EntropyPhaseEvaluator.relative_velocity in all three regimes. The
  tolerance is rtol 2e-4 above a porosity of 0.05; below that it is set by
  Aragog's soft porosity clip.

• The drainage integral is checked on a 60 km martian front entered at
  phi_c = 0.3 (drho = 330 kg/m3, g = 3.711 m/s2), crossed at 70 and
  5 cm/yr (t_res = 85.7 kyr and 1.2 Myr), for five parameter sets:

  grain [m]	eta_melt [Pa s]	mu_s [Pa s]	tau_D / t_res	tau_s / t_res	F_vol at 70 cm/yr	F_vol at 5 cm/yr
  1e-4	10	1e20	8.1653	0.50314	0.27165	0.18134
  1e-3	100	1e22	0.81653	50.314	0.29410	0.22713
  1e-3	100	1e18	0.81653	0.0050314	0.19422	0.10897
  1e-3	1	1e18	0.0081653	0.0050314	0.069535	0.024485
  1e-2	1	1e18	8.1653e-5	0.0050314	0.0094342	0.0025379

  The checks are:

  • The timescales are checked against their closed forms,
    tau_D = L eta_melt / ((1 - phi) drho g (5/7) a^2 phi^4.5) and
    tau_s = mu_s / (drho g L).
  • The drainage equation is separable, so F_vol solves
    t_res = int_F^phi_c max(tau_D, tau_s) dphi / phi. The module matches
    this quadrature, which is independent of the ODE solver, to 1e-6 in all
    ten fractions.
  • In the three cases that stay in the Rumpf-Gupte regime, the module also
    matches the series solution to 1e-5.
  • The matrix-limited case matches phi_c exp(-t_res / tau_s) to 1e-6.

• A dense-melt profile run through the whole front path takes the guard.

• The config validator is tested on test doubles and on a real loaded
  config.

Mutation checks. The new tests fail when any of the following is
introduced:

• releasing the remelted share of the whole mantle instead of the solid, or
  skipping the element reservoirs on release;
• releasing the mean of the solid instead of the ledger's most recent layers,
  popping the oldest layer first, or dropping the fallback for mass the
  ledger does not hold;
• releasing helium twice or not at all, or dropping the noble-gas mole
  updates on restore or on desiccation;
• burying a species whose element total is zero, or raising the closure's
  absolute tolerance to 1e16 kg;
• skipping the desiccation re-check after a remelt, in the wrapper or in the
  main loop;
• setting the totals to the trapped mass on desiccation, emptying the melt
  there, or dropping the mole columns it writes;
• applying the escape floor to the whole total instead of the reachable
  mass, or keeping it on during a remelt, both also in the main loop, or
  ignoring the sign of the release;
• scaling a frozen atmosphere by M_atm instead of the element sum, in the
  unit test and in the main loop;
• emptying the melt on desiccation with trapping on, or booking nothing;
• draining under the surface gravity instead of g_b, a wrong melt density,
  a mobility scaled by 1.1, or not storing the basic-node fields;
• guarding the front on its advance over one step;
• leaving the noble gases out of the buried species;
• emptying the element totals on desiccation whatever the mode, or checking
  the closure by default;
• keeping the trapping columns among the core columns, or not zero-filling
  them;
• dropping (1 - phi) from w_D, taking the faster of the two processes, or
  changing the Rumpf-Gupte exponent from 4.5 to 4.4;
• dropping the withhold, the restore, or both around the chemistry solve;
• removing the supply cap;
• flooring the density contrast in the porosity, or dropping the single-phase
  porosity rule.

Suite results:

• Merged with main up to Add the Aragog core evolution model as core_bc = core_module #980. The helpfile registry, the golden trajectory
  and the migration test keep the trapping columns next to main's accretion
  and core-module columns. The trapping sites are declared in the scanner's
  per-function override table, and the output and config references are
  regenerated.
• Full PR unit tier after the merge with Add the Aragog core evolution model as core_bc = core_module #980, with fwl-aragog 26.10.10:
  5329 passed, none failed on a clean checkout. The golden run passes.
  The coverage figure below was measured at the merge with Turn on Zalmoxis Anderson acceleration by default and require aragog 26.10.08 #971 (5105
  passed).
• Changed-line coverage from the unit tier alone is 96.7% (636 of 658
  executable lines).
• ruff check, ruff format --check, tools/validate_test_structure.sh, the
  four reference generators with --check, and the AGENTS.md checks pass.
  tools/check_test_quality.py --check gives the same per-rule counts as
  origin/main.
• The golden run (tests/integration/test_integration_golden_run.py)
  reproduces every existing column of main's recorded trajectory.
• The porosity rule was checked on the snapshots of an earlier coupled
  Earth-analogue run (1 M_E, 1.5 AU, Aragog with 80 levels, PALEOS). It
  equals the exact volume fraction to 1e-14 at all 19073 mush node samples.

Coupled validation. Earth-analogue runs (1 M_E at 1.5 AU, Aragog with 80
levels, PALEOS, AGNI, CALLIOPE, rfront_loc = 0.4), each a front run and a
trap_mode = "none" control identical otherwise, at 63a808e.
freeze_volatiles is on and the solidification stop is off, so each run goes
on as a frozen mantle under escape until another criterion stops it.

With ZEPHYRUS escape:

	control	front
Time to Phi_global = 0.05	2.01 Myr	1.62 Myr
Surface pressure there	404 bar	276 bar
H in the solid mantle at the end	0	41.2%
O in the solid mantle at the end	0	35.2%
C in the solid mantle at the end	0	1.5%
Frozen rows	133, to 29 Myr	132, to 39 Myr
Escaped while frozen	8.5e19 kg	1.1e20 kg
Stopped on	energy balance	energy balance

• Per-element closure holds to 2e-10 or better on every row of both runs
  (H, C and O). M_atm is 1.000017 to 1.000031 times the element sum in the
  atmosphere.
• Phi_global never rises, so no step remelts.
• 131 front steps had no front while melt remained, all before Phi_global
  fell to 0.46. The crystals still take up D_Z on them: 0.3% of the burial.
• 65 steps drained through the integral (percolation-limited throughout,
  F_tl 0.13 to 0.16), between 2.9e5 and 5.5e5 yr.
• From 5.5e5 yr (Phi_global = 0.30) until the mantle is fully solid, every
  step took the no-drainage bound. 278 of these 279 steps have melt denser
  than the solid inside the front, 115 of them with the front also reaching
  the surface node, and one has the mush split into three layers.
• Those guard steps carry 85% of the buried mass. The dense melt is the
  low-pressure PALEOS behaviour noted under Known limitations, so the trapped
  inventory in this case depends mostly on how dense melt is treated. This is
  the main point for review.
• Burial goes on after the freeze while the last 5% of the mantle
  crystallises, until Phi_global reaches 0 at 2.4 Myr. It moves dissolved
  mass from the melt reservoir into the solid and leaves the frozen atmosphere
  alone. It raises the trapped H from 34.6% to 40.6% of the total. In the
  control the melt reservoir keeps its 2.1e20 kg of H from the freeze on, as
  on main.
• No front step advances more than 0.22 of its own thickness. The exception
  is two steps whose front is a single node of zero thickness, where
  trap_front_courant reads inf and the step takes the guard.

Frozen mantle under strong escape. The same pair with 3e16 kg of N and
the dummy escape at 1e7 kg/s.

	control	front
Time to Phi_global = 0.05	2.00 Myr	1.62 Myr
Surface pressure there	341 bar	227 bar
Frozen rows	118	127
Escaped while frozen	2.7e21 kg	1.8e21 kg
N total at the freeze	1.86e16 kg	1.90e16 kg
N total at the end	1.4e13 kg	6.2e13 kg
Largest closure error	N 3.3e-10, H 1.9e-10	N 3.7e-10, H 1.9e-10

• Escape takes each frozen atmosphere to about 1e-3 bar, where the
  atmosphere drops below mass_thresh. Below that escape stops, as on main,
  just above p_stop = 1e-3 bar, and both runs went on to the 5 Gyr limit.
• The N total is below mass_thresh from 4.4 Myr in the front run and from
  6.1 Myr in the control. Escape draws it down to what the melt and the solid
  hold.
• The front run ends with 3.2e20 kg of H and the control with 2.0e20 kg:
  what trapping buried never reached the atmosphere, so escape could not take
  it.
• Both frozen states were resumed with the energy-balance stop off, once at
  63a808e and once at 2e04fb2, the commit before the fix:
  • the front at 2e04fb2 ran ten steps past the resume, to 7.3 Myr, and
    stopped on the next on Per-element reservoir closure failed for N ... (relative difference 0.0127%), with N at 5% of its value at the freeze
    and C at 7.3e-5;
  • from the same state, the front at 63a808e ran to the 5 Gyr limit with
    every element within 3.7e-10, and ended where the original run did
    (6.25e13 kg of N, 1.1e-3 bar);
  • the control at 2e04fb2, where the check is skipped, ran on to the 5 Gyr
    limit. Its N gap passed 1e-4 at 10.2 Myr, and it ended with 2.1% of its N
    total and 0.3% of its C unaccounted for. At 63a808e the same control
    closes every element to 4e-10.

Test configuration: Habrok cluster (AlmaLinux 9.8, x86_64), Python 3.12
conda environment, with fwl-aragog 26.10.8, fwl-io 26.10.7.1, MORS 26.10.6,
ZEPHYRUS 26.7.24, Zalmoxis 26.9.21 and AGNI at its pin. CALLIOPE, the
outgassing backend, is at 1e3ad73, older than its 26.10.05 pin, so the coupled
runs carry no noble gases. atmodeller is not installed there, so the graphite
path is tested through a stand-in with atmodeller's write convention.

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
• I have updated the docs, as appropriate
• I have added tests for these changes, as appropriate
• I have checked that all dependencies have been updated, as required

[codecov]

Codecov Report

❌ Patch coverage is 93.50850% with 42 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.52%. Comparing base (68281f6) to head (7c4726c).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
src/proteus/outgas/trapping.py	93.73%	9 Missing and 11 partials ⚠️
src/proteus/outgas/compaction.py	90.83%	7 Missing and 5 partials ⚠️
src/proteus/outgas/common.py	78.94%	2 Missing and 2 partials ⚠️
src/proteus/outgas/wrapper.py	95.00%	2 Missing and 1 partial ⚠️
src/proteus/utils/coupler.py	92.00%	1 Missing and 1 partial ⚠️
src/proteus/interior_struct/zalmoxis.py	66.66%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #928      +/-   ##
==========================================
+ Coverage   94.06%   94.52%   +0.45%     
==========================================
  Files         132      135       +3     
  Lines       21681    22284     +603     
  Branches     3644     3745     +101     
==========================================
+ Hits        20395    21063     +668     
+ Misses        888      875      -13     
+ Partials      398      346      -52     

Flag	Coverage Δ
unit-tests	92.85% <93.50%> (+0.05%)	⬆️

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, this is a substantial piece of work, and the bookkeeping closes in a real coupled run with the front engaged. The global front scheme is fine for this PR; per node can follow later. Five things need fixing before this can go in.

Make remelting the reverse of burial. Release is sized from the current melt concentration C_Z, not from what was buried (trapping.py:785-799). Freeze from Phi 1 to 0.5, degas the melt by 100x, remelt to Phi 1: 98 % of the buried water stays in H2O_kg_trapped with no solid mantle left; the kg_liquid <= 0 and melt_mass <= 0 exits release nothing. Size the release from the trapped reservoir (or a burial ledger) and add a test where Phi returns to 1 after C_Z has changed.

trap_mode = 'none' must not change behaviour. run_desiccated now calls keep_only_trapped_mass unconditionally, so on desiccation every H, C, N, S and noble _kg_total becomes 0 where main keeps them. Reproduced on both trees with the same row.

Remove the dt dependence of F_tl. The Courant guard switches to the no-drainage bound when v_f dt / L > 1, and F_tl then rises by about 50 % for a longer step with the same front. t_res = L / v_f does not depend on dt, so the drainage integral stays valid there.

Keep existing runs resumable. A helpfile written by main raises HelpfileSchemaDriftError on this branch, also with trap_mode = 'none'. Give the 32 new columns a default of 0 on read.

Drop the references to the compaction note. Nobody reading this PR can open it. The description, docs/Validation/outgas/compaction.md and test_compaction.py:261-290 should state the equations and expected numbers themselves, and the test should pin the model as shipped, not a variant with the (1 - phi) factor patched out.

Also: noble gases dissolved in the melt are never buried (run_trapping loops over vol_list only), and 21 inline comment blocks are over the 3-line cap.

[IKisvardai]

Thanks Tim, all seven points are addressed in the four commits just pushed. The PR description is updated to match.

1. Remelting is now the reverse of burial. A remelting step returns the share (Phi(t) - Phi(t-1)) / (1 - Phi(t-1)) of every trapped species and element, taking it from the trapped reservoir. Your case is now a test: freeze from Phi 1 to 0.5, degas the melt 100x, then remelt.

   • Remelting to Phi 0.75 returns half.
   • Remelting to Phi 1 returns everything, and _kg_trapped and _kg_solid go to zero.
   • The test also pins the old rule, which would have kept 98.0%.
   • Remelting from a fully solid mantle, where there is no melt to size anything from, releases its share too.
   • Remelting steps now report branch 7.

   The scheme keeps no burial ledger, so the trapped mass is treated as spread evenly through the solid. This is listed as a limitation.

2. trap_mode = 'none' behaves as on main.

   • run_desiccated calls keep_only_trapped_mass only when trapping is on.
   • The per-element closure in assert_mass_conservation now sits behind check_element_closure. It is off by default, and the main loop sets it only when trapping is on.
   • A new test checks that under 'none' desiccation keeps every _kg_total.

3. F_tl no longer depends on dt. The Courant guard is gone; trap_front_courant is still reported. A test runs the same front at the same speed over a step 5x longer (Courant 1.99) and gets the same F_tl and t_res.

4. Runs written by main resume. The 37 trapping columns (trap_* and *_kg_trapped, now including the noble gases) are not core keys. They are filled with 0 on read, with a log line; any other missing column still raises HelpfileSchemaDriftError.

5. No more references to the compaction note. The PR description, docs/Validation/outgas/compaction.md and the test now state the equations and the expected numbers themselves.

   • The test pins the model as shipped, (1 - phi) included.
   • The timescales are checked against their closed forms.
   • The retained fractions are checked against an independent quadrature of the separable drainage equation (t_res = int max(tau_D, tau_s) dphi / phi, to 1e-6), and against its series solution in the Rumpf-Gupte regime (to 1e-5).
   • The retained fractions are pinned at 1e-4.

6. Dissolved noble gases are buried at F_tl C_Z dM_RM. They are withheld from and restored to the chemistry like the other elements, and have _kg_trapped columns of their own.

7. Every inline comment block the branch adds is now at most 3 lines.

The Earth pair was not re-run: none of the changed paths occurs in it. Phi_global never rises, no step has trap_front_courant > 1, neither run carries a noble gas, and neither desiccates.

[IKisvardai]

Thanks Tim! The slow remelt is fixed and tested through the main loop, and I took your suggestion for the closure tolerance. Main is merged up to #947, and the PR description is updated.

Slow remelt after desiccation (162350c). With trapping on, desiccation now changes no total, the same as with trapping off. It empties the atmosphere, keeps each species' and element's trapped mass in the solid, and holds the rest of the total in the melt.

• Mass a remelt releases stays in the melt and the totals until check_desiccation finds it above mass_thresh. That takes two steps in your case.
• With your numbers (H 4e18 kg, O 3.2e19 kg trapped, Phi 0.30 to 0.50), 1000 steps of 2e-4 and 4 steps of 0.05 now end with the same totals, 5/7 of the water still trapped and 1.03e19 kg in the melt.
• The remainder below mass_thresh at the onset of desiccation is now kept in the melt too, rather than deleted.

A second deletion of the same kind (d82096a). Writing the main-loop test turned up another one, in the escape step. Its desiccation floor emptied an element to its locked mass whenever the reachable remainder fell below mass_thresh. Once the flag had cleared, the hydrogen a slow remelt released each step sat below the threshold, and every escape step deleted it again, even at an escape rate of zero. The floor now applies, as on main, when the whole-planet total falls below mass_thresh, and it never goes below the locked mass. The desiccation gate already reads the reachable inventory, so nothing else depended on the old condition.

Main-loop test (d82096a).

• test_a_desiccated_planet_remelting_in_small_steps_keeps_what_it_releases runs Proteus.start on input/dummy.toml. The real dummy interior runs every step, and only Phi_global is set around it.
• The planet is desiccated with 3.6e19 kg of water trapped, then remelts from 0.30 to 0.50 in 8 steps, each releasing less than mass_thresh.
• The totals hold on every row, and the run ends with the trapped 5/7 and the outgassed 2/7 that 4 steps of 0.05 reach.
• The test fails with the old desiccation rule, with the old escape floor, and with the desiccated_after_trapping call removed from the loop.
• It takes about 3 s, or 5 s under coverage.

Closure tolerance (f5e736b). The check no longer uses mass_thresh. It refuses a mismatch once it exceeds both a fixed 1 kg and the solver's relative tolerance times the total. The 6.7e-2 kg of dust N from your Earth run passes; a 2 kg break on a dust total, or a 1e15 kg break on a 5e14 kg total, is refused.

[timlichtenberg]

Thanks Imre, both points are fixed: the slow remelt no longer depends on the step length, and the closure check no longer depends on mass_thresh.

Edit, 4 Oct: #952 and #953 are on main now and take work off this PR, so I rewrote the rest of this comment. Please merge main first; it conflicts in escape/wrapper.py, outgas/wrapper.py, output.md, golden_run.tsv and test_proteus.py.

What main does for you now

• The closure stop on a crystallized planet (3e16 kg of N, freeze_volatiles, trapping on) came from main: the floor zeroed a total whose atmosphere column still held mass. #953 keeps the floor off on a frozen mantle, so there is nothing to fix here anymore. Please run that case again after the merge.
• Desiccation keeps the *_kg_solid columns, sets each total to its solid share and books the removed mass in a new column M_desiccated. check_desiccation leaves that column out of the loss, so a resume gives the same verdict. Trapped mass survives desiccation without extra handling in the trapping code.
• Aragog 26.10.03 (#952) has mobility_function in aragog.eos and writes porosity_b, rho_solid_b, rho_melt_b and g_b at the basic nodes. The transcribed permeability law in compaction.py and the read through _lookup_at_phase_boundary in trapping.py can go.

Still open: the escape floor depends on the trapped mass. On a molten mantle the floor still acts. Two runs that differ only in D_const_N2: with 5e14 kg of N trapped, the floor deletes 1.2e16 kg of reachable N in one step, not booked as escape; with 1.1e16 kg trapped, the total never falls below mass_thresh and the remainder is kept. check_desiccation tests the reachable mass, so the two disagree on what depleted means. Use the reachable mass in both.

The (1 - phi) in darcy_velocity stays, but the reason in the compaction.py docstring is not right. Both forms are the melt velocity relative to the matrix. McKenzie (1984, Eq. 3.7) takes the buoyancy against the mixture, (1 - phi) (rho_s - rho_f) g, and Abe (1995, Eq. 40) takes it against the solid. Please reword that bullet and state the range: McKenzie's solution assumes touching grains, so melt fractions below about 0.2, and the front enters at rfront_loc = 0.4 by default. FormingWorlds/aragog#67 tracks the same factor for Aragog.

Two comments still state the opposite of the code (escape/wrapper.py 607, outgas/wrapper.py 335).

I can send the configs of the runs.

[IKisvardai]

Thanks Tim. Main is merged (up to #954), and all your points are in. The PR description is updated.

Merge with #952 and #953 (18fb224).

• Desiccation now follows Keep the escape floor off a frozen mantle and book desiccation #953: the atmosphere is emptied and booked in M_desiccated, and the solid share is kept, trapped mass included. The branch's hold_in_mantle is gone.
• One deliberate difference: with trapping on, run_desiccated also keeps the melt. Without that, a slow remelt's release would be removed and booked every step, and what the planet keeps would depend on step length again, the point of your last review. With trapping off it is exactly Keep the escape floor off a frozen mantle and book desiccation #953.
• _solve_chemistry passes the new initial flag on.
• The escape module now imports the trapping helpers inside the functions that use them, to break an import cycle with outgas.wrapper.

Crystallised N case, after the merge. I reproduced the path with the dummy modules, a trapping run and a trap_mode = "none" control:

• 2.9e16 kg of N, 3.9e17 kg of H, freeze_volatiles on;
• the mantle frozen at Phi 0.002 before escape (bulk reservoir) starts.

On 37 frozen steps the N total is below mass_thresh while the atmosphere still holds 6.6e13 to 6.5e15 kg of it. Both runs complete, every element closes to 4e-15 relative on every row, and the two agree. If you send the configs I'll run your case as it is.

Escape floor on the reachable mass (3179257). The floor now tests the total less what trapping locks, which is what check_desiccation reads, and empties it down to the locked mass.

• A test runs your two cases, 5e14 and 1.1e16 kg of N trapped. Both now empty the reachable remainder; the whole-total floor kept it in one and not the other.
• On its own this would bring back the step-length dependence from your last review. Once a desiccated planet's flag clears, a slow remelt's released H sits below mass_thresh and would be floored every step.
• So run_escape turns the floor off on a step on which a remelt returned buried mass (trap_branch == 7, trap_kg_step < 0), the same way Keep the escape floor off a frozen mantle and book desiccation #953 turns it off on a frozen mantle.
• The main-loop test (now 16 steps of 0.0125) fails without that switch, and fails with the whole-total floor.

Aragog 26.10.03 (727b4d6, 4592d2f).

• compaction.mobility_function now calls aragog.eos.mobility_function, keeping only the grain-size check, and the transcription is gone.
• The aragog interior stores phi_basic, rho_basic, rho_solid_b, rho_melt_b and g_b, and the front is located and drained on the basic nodes. The reads of _lookup_at_phase_boundary, phase_staggered._g and the structure profile are gone. The lowest basic node is the CMB, so the special case for a front resting on it went too.
• I did not use porosity_b for the drainage. Its docstring says it reads about 0 where rho_s <= rho_l, and the Earth front had dense melt on 257 steps, which carried 81% of the burial; it is also not the melt fraction at a single-phase node. The porosity stays the sign-keeping lever rule, now on Aragog's densities. A new test pins it to Aragog's _porosity_from to 2e-5 where the melt is lighter, the soft clip being 2.5e-7/phi.
• One consequence: the front geometry changes in detail, so the Earth-pair numbers in the description need a re-run before they are quoted. The description now says this.

(1 - phi) (727b4d6). The docstring bullet now says what you wrote: both forms are the velocity relative to the matrix, with McKenzie (1984, Eq. 3.7) taking the buoyancy against the mixture and Abe (1995, Eq. 40) against the solid. It states that McKenzie's solution holds below about 0.2 while the front enters at rfront_loc = 0.4, and points to FormingWorlds/aragog#67. I took the equation numbers from your comment.

Comments. Both are rewritten. escape/wrapper.py now describes the reachable floor. The desiccated_after_trapping docstring now says the release is kept in the melt until the flag clears.

[timlichtenberg]

Thanks Imre, this round closes most of it: the slow remelt after desiccation, the floor on the reachable mass and the drainage on Aragog's exports all hold, also in runs. I hope the point below is the last change before this can go in.

A frozen-mantle trapping run with escape still stops, now on the closure check. With the R3 setup from my last review the run gets past the floor, but raises in the per-element closure check at about 2.7e4 yr.

The mismatch behind it is tiny and already on main: escape debits each element total by its share of the escaping mass, while run_crystallized scales the atmosphere by (M_atm - esc)/M_atm, and M_atm is 1.0000177 times the sum of the element masses in the atmosphere. That leaves about 1.8e-5 of each step's N loss unaccounted, which is harmless in itself, but as escape draws the N total down the relative gap keeps growing and crosses the 1e-4 limit, so any long enough run of this kind ends in the RuntimeError. With trapping off the same drift runs on, because the check is skipped.

Please make the two debits remove the same mass, for example by debiting each element's atmosphere by what escape took from its total, or let the check allow for that drift, and add a frozen-mantle escape case over enough steps to the main-loop test.

For the next rounds, and future PRs that touch the volatile bookkeeping: please run a few coupled test cases to the end before pushing, including a frozen mantle with escape over a long stretch. Most of the points in the last rounds only showed up after many steps, and a full run catches them much earlier than the unit tests can.

On the reply: the main-loop test passes with the whole-total floor in my run; only the unit test fails without the switch.

[IKisvardai]

Thanks Tim. The frozen-mantle drift is fixed. Two coupled pairs ran to the end before the first push; a third was still running then. I'll do the same before every push that touches the volatile bookkeeping. The PR description is updated.

Frozen-mantle escape (42be63c). I took your first option: the two debits now remove the same mass.

• run_crystallized now takes the retained fraction over the element masses in the atmosphere, the same sum calc_new_elements shares the loss over, instead of M_atm. So each element's atmosphere loses exactly what escape took from its total.
• The species columns and pressures scale by the same fraction, which keeps the composition and the ratio of M_atm to the element sum. Keeping the species at the M_atm-based fraction would instead have held their absolute difference fixed while both deplete, leaving a phantom species atmosphere once the elements are gone.
• A row that records no element atmosphere keeps the species sum (only some unit-test rows do this).
• The fix acts with trapping off as well. The golden run is unchanged, since it never escapes from a frozen mantle.

Tests:

• A unit test runs 60 alternating escape and frozen steps on an H2O-N2 atmosphere at your 1.0000177 ratio. Every element closes to 1e-10 while the N total falls to what the solid holds; under the old scaling the N gap would be 2.1e-4 of the total.
• The main-loop test you asked for (test_a_frozen_mantle_escapes_over_many_steps_with_every_element_closing) freezes the dummy mantle with freeze_volatiles and trapping on, then escapes a quarter of the atmosphere a step for 20 steps. The dummy chemistry is wrapped to put the species masses at the CALLIOPE ratio, since its own split is exact.
• With the old scaling that loop stops on the closure check, Per-element reservoir closure failed for H ... (relative difference 0.0107%), the same failure as your R3 run. With the fix every element closes to 1e-10 on every row.

Correction to my last reply. You're right: the main-loop test passed with the whole-total floor, and my claim that it failed was wrong. My mutation harness copied only src/ into a scratch tree. get_proteus_dir then found no pyproject.toml beside it, so every main-loop mutant errored in Proteus.start whatever was mutated, and counted as a kill. The same affected the main-loop claims in my third-round reply. I've fixed the harness and re-run them all:

• the remelt floor switch, the kept melt and the desiccation re-check are genuinely caught by the main-loop test;
• the whole-total floor was not. I've added 1.5e18 kg of trapped N2 to that test (63a808e). Its released 2/7 stays below mass_thresh, is kept while the remelt lasts, and is floored to the locked 5/7 once it stops. That now fails under the whole-total floor.

Coupled runs. All at 63a808e on the Earth analogue from the description, each with a trap_mode = "none" control. freeze_volatiles is on and the solidification stop is off, so every run goes on frozen under escape until another criterion stops it. The numbers are in the description.

• With ZEPHYRUS escape, both runs ended on energy balance, the front at 39 Myr after 132 frozen rows and the control at 29 Myr after 133. At 1.5 AU escape takes little from the frozen atmospheres (1.1e20 and 8.5e19 kg). Every element closes to 2e-10 on every row. The description's Earth-pair numbers now come from this pair, with the front on the basic nodes. At Phi_global = 0.05 the front run has 34.6% of the H, 29.4% of the O and 1.4% of the C trapped at 1.62 Myr, against 35.1%, 29.8% and 1.4% at 1.59 Myr on the staggered nodes.
• For a long stretch of strong escape, the same pair with 3e16 kg of N and the dummy escape at 1e7 kg/s. Escape took both frozen atmospheres from a few hundred bar to about 1e-3 bar, over more than 100 frozen rows each, and both ran to the 5 Gyr limit. N fell from 1.9e16 kg to 6.2e13 kg in the front run and 1.4e13 kg in the control, below mass_thresh most of the way. Every element closes to 4e-10 on every row of both.
• I resumed those frozen states at 2e04fb2, the commit before the fix, with the energy-balance stop off. The front ran ten steps past the resume, to 7.3 Myr, and stopped on the next on Per-element reservoir closure failed for N ... (relative difference 0.0127%), with N at 5% of its value at the freeze. That is your failure, in a coupled run. The control at 2e04fb2 ran to 5 Gyr with 2.1% of its N total and 0.3% of its C unaccounted for, and nothing flagged it. At 63a808e the same resumes close every element to 4e-10.
• One thing these runs show that the earlier pair, stopped at Phi_global = 0.05, did not: burial goes on after the freeze while the last melt crystallises, until Phi_global reaches 0. It moves melt into the solid and leaves the frozen atmosphere alone. In the ZEPHYRUS pair it takes the trapped H from 34.6% to 40.6% of the total. It is now listed under modelling choices.
• A third pair with rfront_loc = 0.9 is still cooling: Phi_global is 0.22 at 5.3 Myr after 8 h of wall time. Every element closes to 2e-10 on its 1400 rows so far, and Phi_global has never risen. Both runs of this pair, the control included, show episodes where the surface warms by up to about 400 K within about 10 kyr while Phi_global keeps falling. The other two pairs show none, so the episodes come with rfront_loc = 0.9 in the interior, not from trapping. I'll post the pair's end when it gets there.

Environment for these runs. The cluster moved to AlmaLinux 9.8 this week.

• MORS, ZEPHYRUS, fwl-aragog (26.10.08) and fwl-io (26.10.7.1) are at main's pins.
• CALLIOPE is at 1e3ad73, older than the 26.10.05 pin, because I'm keeping that checkout unchanged. That version has no noble-gas support, so I couldn't run your He and Ar case. The two Seed the CALLIOPE cold start so identical runs agree #958 cold-start smoke tests fail against it.
• The datasets for the new fwl-io layout were fetched for these configs rather than relocated.

Happy to run your R3 config directly if you send it.

[timlichtenberg]

Thanks Imre. The frozen-mantle drift is closed, the mole columns follow the masses, the three unused columns are gone and main is merged. Approved; please merge only once all checks on this head are green.