feat(antarctica): ISMIP7 core-experiment pipeline with n=3 rheology and hist-branched CTRL

[hoffmaao]

Intent

Correct the rationale for the hist-branch CTRL change (comment + operator-warning text only; no behavior change). Context: an earlier commit (01696a7, already on PR #4) justified branching the CTRL from the historical endpoint by claiming the CTRL/projection apparent-MB (a_ref) mismatch of -95 vs -777 Gt/yr was a constant net sink that leaked a spurious ~160 mm SLE trend into proj-CTRL, and that removing it dropped us into the ISMIP6 envelope. That reasoning was WRONG: a_ref is a t=0 balancing correction that zeroes the initial thickness tendency, not a free net sink, so multiplying 681 Gt/yr by 85 years estimates nothing physical. An isolation test at fixed n=4 rheology (same ssp126 projection, only the control baseline swapped) shows the cold-start CTRL (its own drift +8 mm) gives proj-CTRL=+132 mm while the hist-branch CTRL (drift +37 mm) gives +161 mm - i.e. the CORRECT same-state control makes the signal LARGER, because the projection's historical-endpoint relaxation only cancels when the control shares it. The hist-branch change is KEPT (it is the protocol-correct same-state differencing the ISMIP6 ctrl_proj convention requires), but the ISMIP6 overshoot is now understood to be a real forced-response bias (melt/dynamics), not a differencing artifact; n=3 rheology roughly halves it (ssp126 +85 mm, still above the [-14,+50] envelope). This commit only rewrites the code comment and the cold-start warning message to state the correct rationale. This lands on antarctica-n3 / PR #4, which stays OPEN (Shapero merges to main later).

What Changed

• Builds the ISMIP7 Antarctic simulation framework: icepack2_tools gains forcing readers (ISMIP7 atmosphere/ocean, RACMO baseline), regridding, dual friction, and a thermal model; antarctica/scripts adds the shared experiment runner, control/projection drivers for all core experiments (historical, ssp126/370/585 x 2 ESMs, CTRL, OCX), conservative thickness transport with mass-budget audits, rank-robust warm restarts, a solver rescue ladder with per-step dt-subcycling, and gated when-ready launchers with preflight checks.
• Adds the n=3 rheology track: a physical-prior inversion (thermo-derived fluidity plus a Whittle-Matern prior) with a converged, validated MAP, n-tagged MAP filenames threaded through preflight, launchers, and the forward drivers, and loud failure on non-converged inversion solves. The CTRL is now branched from the historical endpoint (protocol-correct same-state differencing), with the rationale corrected in code comments and docs: the ISMIP6 envelope overshoot is a real forced-response bias, not an a_ref differencing artifact, and n=3 roughly halves it.
• Records results and tooling for the matrix: 9/11 core experiments complete at 32 km with per-core reports, MATRIX_STATUS.md, an ISMIP6 ensemble comparison script with scenario-aware pool selection, an ISMIP6-track audit script, per-basin melt calibration, and a reworked Globus forcing mirror for the reorganized GHub share. The no-mistakes review pass surfaced 30 findings; the two documentation ones (stale hist-branch rationale in README/MATRIX_STATUS and the cores 9/10 cold-start CTRL provenance note) were fixed on-branch, and the remaining code-level findings are tracked as open items.

Risk Assessment

✅ Low: The new delta is a four-line docs-only addition to MATRIX_STATUS.md that implements the approved provenance-note fix verbatim in spirit - factual, no retracted claim reintroduced, no code or behavior change - leaving nothing new to flag.

Testing

Proved the change is behavior-neutral via normalized-AST comparison of run.py before/after, then exercised the real CTRL driver in the Firedrake environment through both branch paths to capture the actual operator-facing messages, confirming the corrected same-state-differencing rationale replaces the retracted a_ref-sink claim in code, warning text, and docs; all checks passed and the worktree was left clean.

Evidence: CTRL operator messages (live transcript, both branch paths)

CASE A: no historical endpoint present -> cold-start WARNING path WARNING: no historical endpoint hist_cesm2_waccm_evharness_2500_final.h5; cold-starting from the inversion. The CTRL will start from a DIFFERENT geometry than the hist-branched projections, so projection-minus-CTRL will not cleanly isolate the forced response. Run the historical first, or pass --restart. CASE B: historical endpoint present -> hist-branch message path Branching CTRL from the historical endpoint (same initial state as the projections; shared drift cancels in proj-CTRL): hist_cesm2_waccm_evharness_2500_final.h5

========================================================================
CASE A: no historical endpoint present -> cold-start WARNING path
========================================================================
  WARNING: no historical endpoint hist_cesm2_waccm_evharness_2500_final.h5; cold-starting from the inversion. The CTRL will start from a DIFFERENT geometry than the hist-branched projections, so projection-minus-CTRL will not cleanly isolate the forced response. Run the historical first, or pass --restart.

========================================================================
CASE B: historical endpoint present -> hist-branch message path
========================================================================
  Branching CTRL from the historical endpoint (same initial state as the projections; shared drift cancels in proj-CTRL): hist_cesm2_waccm_evharness_2500_final.h5

harness done; fake endpoint removed: True

Evidence: Harness used to drive run.py:main() through the branch decision

"""Drive antarctica/scripts/control/run.py:main() through the CTRL-branching
decision for real, stopping just after the operator message (before any mesh
or solver work) by stubbing setup_model."""
import os, sys

WORKTREE = os.environ["WORKTREE"]
sys.path.insert(0, os.path.join(WORKTREE, "antarctica", "scripts", "control"))

import run  # noqa: E402  (imports firedrake + the real module under test)


class StopAfterBranch(Exception):
    pass


def fake_setup_model(*a, **k):
    raise StopAfterBranch


run.setup_model = fake_setup_model

esm_tag = run.ESM.lower().replace("-", "_")
tag = "evharness"
hist = os.path.join(
    run.RESULTS_DIR, f"hist_{esm_tag}_{tag}_{run.lc}_final.h5"
)

print("=" * 72)
print("CASE A: no historical endpoint present -> cold-start WARNING path")
print("=" * 72)
assert not os.path.exists(hist), f"unexpected pre-existing file: {hist}"
sys.argv = ["run.py", "--tag", tag]
try:
    run.main()
except StopAfterBranch:
    pass

print()
print("=" * 72)
print("CASE B: historical endpoint present -> hist-branch message path")
print("=" * 72)
os.makedirs(run.RESULTS_DIR, exist_ok=True)
with open(hist, "wb") as f:
    f.write(b"placeholder")
try:
    sys.argv = ["run.py", "--tag", tag]
    try:
        run.main()
    except StopAfterBranch:
        pass
finally:
    os.remove(hist)
print()
print("harness done; fake endpoint removed:", not os.path.exists(hist))

Pipeline

Updates from git push no-mistakes

✅ **intent** - passed

✅ No issues found.

⏭️ **Rebase** - skipped

• ⚠️ antarctica/scripts/control/run.py - merge conflict rebasing onto origin/main
• ⚠️ antarctica/scripts/simulation.py - merge conflict rebasing onto origin/main

🔧 **Review** - 30 issues found → auto-fixed (2) ✅

• 🚨 icepack2_tools/forcing.py:755 - Flotation mask density mismatch: both melt callbacks compute haf = s - (b + (_RHO_WATER/_RHO_ICE)max(-b,0)) with freshwater _RHO_WATER=1000 while the model surface s is built with seawater flotation (rho_water=1024 per icepack2.constants and simulation.py), so genuinely floating ice with h between ~0.866depth and ~1.117*depth is classified grounded and receives zero ocean melt - exactly the near-grounding-line band where melt matters most. Same expression at forcing.py:755 (CTRL climatology callback) and forcing.py:824-826 (projection callback). Ask-user because fixing it changes scientific results and the per-basin K calibration may have partially compensated.
• 🚨 antarctica/scripts/run_control.py:303 - run_control.py's VAF and total-mass diagnostics multiply by rho_I imported from icepack2.constants (ice_density in icepack2's MPa-m-yr unit system, not 917 kg/m^3), so the reported mm SLE and Gt values are off by many orders of magnitude; simulation.py fixed the identical diagnostics with _RHO_I_SI=917.0 (lines 1299-1300) but this script was not updated.
• 🚨 antarctica/scripts/simulation.py:357 - u_c is set from a rank-LOCAL mean of observed speed (Function(Q).interpolate(...).dat.data_ro.mean() with no comm reduction, unlike the mesh.comm.allreduce used elsewhere in the file), so under mpiexec each rank builds a different sliding scale in K_base (line 391) and K_linear (line 406) and the assembled momentum operator is partition-dependent - results silently change with rank count, defeating the file's own rank-count-robust restart design. run_control.py:139 has the same rank-local mean.
• 🚨 antarctica/scripts/inversion_icepack2.py:549 - The inversion's solve-failure sentinel returns last_good_obj[0]*10, but last_good_obj stores the misfit-only J_val (line 556) while successful evaluations return total = J_val + reg_theta + reg_phi; with GAMMA defaults of 1e4 the regularization terms can dominate late in optimization, so a mid-run forward/adjoint failure can return an 'objective' lower than the current total with a zero gradient, which L-BFGS-B accepts and then terminates on projected-gradient 0, silently writing the MAP from the failed control point; the sentinel should store the last good TOTAL (both the forward guard at ~549 and the adjoint guard at ~574).
• ⚠️ icepack2_tools/forcing.py:37 - smb_kgm2s_to_myr multiplies by a spurious _RHO_WATER/_RHO_ICE factor (kg/m^2/s * sec/yr / rho_ice already IS m/yr ice equivalent), inflating every ISMIP7 acabf SMB anomaly by ~9%, inconsistent with the RACMO loader which divides only by rho_ice; already tracked as open item 1 in MATRIX_STATUS.md as 'a real pending decision', so it needs the author's call.
• ⚠️ antarctica/scripts/run_control.py:1 - run_control.py is a stale duplicate of the simulation.py pipeline that is scientifically inconsistent with the branch: it loads the untagged (legacy n=4) MAP but applies n=3 Glen, uses the non-conservative legacy transport with mass-injecting h_clamp=10 L2 projection, truncates nsteps with int(), and writes non-restartable checkpoints without t_yr; it is superseded by control/run.py and should be deleted or reduced to a shim before someone runs it as 'the control'.
• ⚠️ antarctica/scripts/simulation.py:1004 - Checkpoints record friction and lc attrs (and restarts hard-fail on friction/a_ref mismatch) but not n_flow/a4_factor/m_slide, so resuming a checkpoint under a different ISMIP7_N_FLOW or ISMIP7_A4_FACTOR (e.g. resuming an n=4 antarctica-branch run on this n=3 branch, whose defaults differ) runs cleanly with silently wrong rheology - the same failure class the friction guard at line 284 was added to prevent.
• ⚠️ antarctica/scripts/inversion_icepack2.py:658 - If the last scipy evaluation ended in a caught ConvergenceError, n_flow/m_slide are left at mid-continuation values (e.g. n=1.5), and the final slvr.solve() can converge at the wrong exponent without raising, so the 'Final misfit (masked)' and the velocity field written into the MAP checkpoint use the wrong rheology; reassign n_flow/m_slide to their target values before the final solve.
• ⚠️ antarctica/scripts/launch_rc_ctrl_when_ready.sh:99 - Presence-vs-value env-var mismatch: launching with ISMIP7_FIXED_FRONT=0 (intending fixed-front OFF) still runs with fixed front because the script only skips its own export while the inherited variable remains set and simulation.py:828 tests presence (is not None); the non-1 branch needs unset ISMIP7_FIXED_FRONT. launch_rc_inversion_when_ready.sh:89 has the same bug with ISMIP7_NO_CALVING_TERMINUS=0 silently dropping the calving-terminus term (inversion_icepack2.py:155 tests presence).
• ⚠️ antarctica/scripts/experiment.py:121 - The pre-run atmosphere check only verifies that SOME acabf years exist, not that they cover [t_start, t_end], and ISMIP7Atmosphere.get_field silently returns zeros for a missing year; in full-field fallback mode (smb_anomaly=False, taken whenever RACMO is absent) any coverage gap silently applies zero total SMB for the missing years, corrupting the mass balance with no warning. Ask-user because a zero-anomaly extension is partly intentional (the tracked 2014-to-2015 handoff gap relies on it), so the right fix is a coverage check against the run window rather than removing the zero fallback.
• ⚠️ antarctica/scripts/experiment.py:113 - Projections have no auto-resume path (ISMIP7_AUTO_RESUME is only wired into control/run.py): rerunning a crashed multi-week projection without a manual ISMIP7_RESTART restarts from hist_*_final.h5 (t=2014), and run_simulation's CSV truncation (simulation.py:1041-1053) then deletes every timeseries row after 2014.5, silently discarding all of the crashed run's progress; wiring latest_checkpoint into experiment.py would mirror the control driver.
• ⚠️ icepack2_tools/regrid.py:56 - ISMIP7Regridder computes its query-point bounding box from rank-local coordinates.dat.data_ro, but Firedrake VertexOnlyMesh defaults to redundant=True (points from rank 0 only), so under MPI the output grid is silently cropped to rank 0's partition bbox and the per-rank _query_shape/crop.reshape mismatches on other ranks - parallel regridded NetCDF output is wrong or crashes.
• ⚠️ icepack2_tools/forcing.py:288 - ISMIP7Atmosphere.get_field silently returns zeros when a year's file is missing (_load_year returns None for any missing or misnamed path), and get_thermal_forcing has the same silent-zero fallback for a missing variable directory, so a renamed file or a run past the last available year quietly zeroes the forcing mid-run instead of failing; at minimum this should warn loudly once per missing year.
• ⚠️ antarctica/scripts/core_report.py:56 - core_report.py always passes its --exps default "exp01..exp05" to compare_ismip6.py, so the scenario-aware auto-pool selection added in commit 01696a7 can never engage via core_report: an ssp126 report is graded against a mostly-RCP8.5 ensemble pool (only exp03 is RCP2.6) and the tracked 'outside envelope' verdicts (core05/06/08 reports) use the wrong comparison set; the default should be None so auto-selection runs, but this changes graded verdicts in tracked reports.
• ⚠️ antarctica/scripts/core_report.py:87 - core_report.py conflates tool failure with a scientific result: any nonzero exit from compare_ismip6.py is written to the tracked report as 'outside envelope', but exit code 2 means usage/data error (missing ~/data/ismip6 archive, no year overlap), and audit rc==2 is similarly recorded as 'OFF TRACK' when the CSV merely lacks budget columns - a broken comparison gets permanently recorded as an out-of-envelope finding.
• ⚠️ antarctica/scripts/download_forcing.py:313 - download_forcing.py skip-if-exists uses only local_path.exists() with no size/checksum comparison, so a file truncated by an interrupted transfer prints '[exists]', is excluded from the TransferData list, and the sync_level="checksum" repair never runs - the corrupt file silently feeds the pipeline indefinitely.
• ⚠️ antarctica/scripts/download_forcing.py:406 - download_forcing.py scenario-group completeness is a bare file-count test (n_local < n_remote), so truncated local files or stale extras make the count match, the group is marked complete and skipped, and the checksum-sync transfer that would repair corruption is never submitted.
• ⚠️ antarctica/scripts/download_forcing.py:394 - download_forcing.py treats a transient Globus auth/network error the same as absent data: list_remote_files returning [] makes _pick_version return None, reported as 'not on the share, skipped', and the run exits 0 - a whole (ESM, scenario) forcing set can be silently omitted from the mirror; a listing failure should be distinguished from an empty listing and exit nonzero.
• ⚠️ antarctica/README.md:268 - Stale rationale: the README still describes the cold-start CTRL warning with the retracted claim ('its apparent-MB baseline will not match hist-branched projections (a spurious trend in proj-CTRL)'), contradicting the corrected rationale in commit 8f873a7 (a_ref is a t=0 balancing correction, not a net sink; the overshoot is a real forced-response bias) and no longer matching the actual warning text in control/run.py:208-214.
• ⚠️ antarctica/reports/MATRIX_STATUS.md:67 - MATRIX_STATUS.md open item 1 (lines 62-67) still asserts the retracted causal claim as fact ('The matrix sitting above the ISMIP6 envelope ... turned out to be the CTRL cold-starting ... with a mismatched frozen apparent-MB baseline' and that cores 9/10 'carry the spurious trend'); per the corrected rationale the overshoot is a real forced-response bias and the hist-branch control makes proj-CTRL larger (+161 vs +132 mm), not smaller, so this status doc now misstates the project's own conclusion.
• ℹ️ antarctica/scripts/calibrate_melt.py:259 - calibrate_melt.py reduces with rank-local numpy sums (node_kgyr[mask].sum(), lumped.sum()) despite the Firedrake/PETSc setup working under MPI, so running it under mpiexec silently produces wrong per-basin K and K* with normal-looking output; an nprocs>1 guard would prevent silent misuse.
• ℹ️ antarctica/scripts/launch_ctrl_when_ready.sh:1 - The three when-ready launchers are ~85% line-for-line identical (noclobber lock protocol, traps, resource readers, gate loop, NTAG awk, logging) across ~315 lines and have already drifted (boundary_ids preflight in 2 of 3, MIN_FREE_CORES sanity warning in 1 of 3); consolidating the gate/lock machinery into one sourced library or a single parametrized launcher would make future fixes to the subtle lock/launch logic land in one place.
• ℹ️ antarctica/scripts/preflight.py:32 - The n-tag MAP-filename rule now has at least 6 implementations that must stay bit-identical (inversion_icepack2.py:94, preflight.py:32, simulation.py:65, plus awk one-liners in all three launchers); a divergence silently makes a launcher gate on, or a forward load, a differently-named MAP - at minimum the three Python copies should import one shared helper.
• ℹ️ antarctica/scripts/launch_rc_inversion_when_ready.sh:72 - Unlike both CTRL launchers, the inversion gate never pre-checks boundary_ids.json (required at inversion_icepack2.py:152), so the gate can wait hours for resources, fire, die in seconds on FileNotFoundError, and retain the lock as an 'already launched' marker that masks the failure.
• ℹ️ antarctica/scripts/simulation.py:790 - Because t_restart unconditionally overrides the caller's t_start, hist-branched projections/CTRL integrate 2014->end (not the nominal 2015->end): year 2014 is replayed under the scenario forcing with a typically absent (zero-filled) 2014 anomaly - consistent between CTRL and projections so it cancels in differences, but it shifts absolute timeseries and duplicates a year already simulated in the historical (related to tracked open item 2 in MATRIX_STATUS.md).
• ℹ️ antarctica/scripts/projections/ocx.py:45 - OCX always cold-starts (setup_model() with no restart argument and no --restart/auto-resume handling), so a reboot of the 1990-2025 run at dt=0.1 silently starts over and rewrites the timeseries CSV from scratch even though valid self-contained checkpoints exist on disk.
• ℹ️ icepack2_tools/coupled.py:139 - coupled.py forcing_callback multiplies melt_rate_ice_equivalent() by 31556926 assuming the plume model returns m/s and copies it dof-for-dof into ctx["ocean_melt"] assuming identical dof ordering, but no plume module exists in the repo to verify either assumption - a m/yr-returning plume would be inflated ~3.2e7x and a submesh would scramble the copy.
• ℹ️ antarctica/scripts/compare_ismip6.py:169 - compare_ismip6.py does not re-baseline 'ours' after clipping to the ensemble overlap window, so a run starting well before ~2015 carries its accumulated pre-window signal as a constant offset against ensemble members baselined at their own t0; harmless for the standard 2015-start runs but wrong for earlier-start inputs.
• ℹ️ .gitignore:3 - The bare 'data' pattern in .gitignore ignores any file or directory named 'data' at any depth, not just the repo-root data dir; '/data/' would scope it to the root.
• ℹ️ antarctica/scripts/control/run.py:179 - Intent conformance verified: commit 8f873a7 touches only the comment block (control/run.py:179-196) and the two PETSc.Sys.Print operator messages; the hist-branch behavior, cold-start fallback, and restart validation are unchanged, and the corrected rationale (a_ref as t=0 balancing correction, +132 vs +161 mm isolation-test numbers, real forced-response bias, n=3 halving) is stated exactly as the intent requires; the retracted '681 Gt/yr / ~160 mm' arithmetic survives only inside this comment where it is explicitly quoted as retracted.

🔧 Fix: Align README and MATRIX_STATUS with corrected hist-branch rationale
2 infos still open:

• ℹ️ antarctica/reports/MATRIX_STATUS.md:63 - The MATRIX_STATUS rewrite correctly retracts the wrong a_ref attribution, but it also dropped the still-true provenance note that cores 9/10 predate the hist-branch fix (their tracked CTRL records are cold-start controls, i.e. not the protocol-correct hist-branched control under the corrected rationale), and neither core09 nor core10 report records this itself - core09 even states 'This is the CTRL that the projections difference against'. Consider re-adding one sentence noting the archived core09/10 CTRLs are cold-start (mismatched-baseline) controls; ask-user since it is wording in the tracked scientific record.
• ℹ️ antarctica/reports/MATRIX_STATUS.md:63 - Fix conformance verified: commit b6a3ccb changes only antarctica/README.md (cold-start warning description now matches the actual control/run.py:208-214 text - 'DIFFERENT geometry ... will not cleanly isolate the forced response' - and the branching sentence now attributes cancellation to shared relaxation drift plus the identical frozen a_ref) and MATRIX_STATUS.md open item 1 (restated as a real forced-response bias with the +161 vs +132 mm isolation test, a_ref as t=0 balancing correction, explicit retraction of the '~160 mm spurious a_ref trend' claim, and the n=3 halving to +85 mm vs the [-14,+50] envelope); no code or behavior was touched and no other files changed, matching both the user's fix instructions and the stated intent.

🔧 Fix: Restore cores 9/10 cold-start CTRL provenance note
✅ Re-checked - no issues remain.

✅ **Test** - passed

✅ No issues found.

• AST comparison of antarctica/scripts/control/run.py at 8f873a7^ vs 9b0e36b with string literals normalized - identical, proving the change is comments/messages only
• git show --stat b6a3ccb 9b0e36b - confirmed the follow-up commits touch only markdown (README, MATRIX_STATUS)
• Ran run.py:main() end-to-end through the CTRL branch decision in the venv-firedrake-2026 environment (setup_model stubbed) for both paths: no historical endpoint (corrected cold-start WARNING printed) and fake historical endpoint present (corrected hist-branch message printed, hist-branch behavior kept)
• Repo-wide grep for retracted-rationale text (spurious, apparent-MB baseline, 681, ~160 mm) - only explicit retractions and one unrelated numerics comment remain
• Cross-checked intent numbers (+8/+37 mm drift, +132/+161 mm proj-CTRL, n=3 +85 mm, [-14,+50] envelope) against run.py comment, README, and MATRIX_STATUS
• git status after harness run - worktree clean, fake endpoint h5 removed

✅ **Document** - passed

✅ No issues found.

✅ **Lint** - passed

✅ No issues found.

✅ **Push** - passed

✅ No issues found.

[hoffmaao]

Moving staging to the fork (hoffmaao/ismip7) while this work matures; will reopen or file a fresh PR when it is ready for upstream review.