v2.7 Multi-Layer Box Model Scoping¶
Status: Scoping + minimum-viable experiment (this PR).
Parent finding: v2.6 GEOTRACES hybrid recovered alpfe + Smallgrow + diatomgraz to Cal-grade across seeds; stuck pair shifted to (scav_rat, R_PICPOC) with scav_rat at 0.92 ± 0.01 off Carroll across all GEOTRACES weights, identified as HNLC-region surface-only structural limit.
Hypothesis: A 2-layer (surface + subsurface) box exposes the vertical DFe gradient that scav_rat controls, breaking the surface-only degeneracy.
0. Where this work fits the broader differentiable-emulator landscape¶
DarwinDiff sits in the same architectural family as the AI weather/climate models recently surveyed in AIMIP Phase 1 (Henn et al. 2026, arXiv:2605.06944): autograd-through-physics in our case, alongside autoregressive full-state emulators (ACE2, NeuralGCM, ArchesWeather), generative-diffusion samplers (cBottle), and hybrid physics/AI cores (NeuralGCM-HRD, DLESyM). AIMIP's headline observation is directly relevant to our framing: "AIWCMs typically produce smaller biases than [the physics-based reference] CM4, despite being forced by approximately the same SSTs." The same logic applies to DarwinDiff vs ECCO-Darwin: a differentiable counterpart to a physics-based calibration is allowed to match-or-exceed on the same observations, not just emulate.
Two further AIMIP findings inform v2.7's honest scope:
- Out-of-sample is the hard test. AIMIP reports "major disagreement in the models' response to the very much out-of-sample perturbed SST experiments," paralleling our HNLC-Eq-Pacific vs cross-basin distinction. Our v2.7 result is, intentionally, the in-sample case for GEOTRACES coverage; the structural-limit claim on
scav_ratis explicitly contingent on this being a local AOI fit and remains testable only at cluster scale. - Architectural diversity is itself a finding. AIMIP submissions span "a range of frameworks (autoregressive full emulation and hybrid physics/AI models as well as conditional diffusion), and AI architectures (UNets, Fourier operators, vision transformers, etc.)." DarwinDiff's autograd-Carroll-6 / DINN-per-cell stack is one valid point in that landscape. The point of these intercomparisons is shared evaluation protocol, not architectural convergence.
These citations are background framing for the eventual manuscript intro; they don't change any v2.7 result. Reference: AIMIP authors include Brian Henn and Christopher Bretherton at Ai2, where the OlmoEarth workshop was hosted.
1. Why a multi-layer model is the realistic next lift¶
v2.6 fits a single surface layer (h_mld = 50 m). The iron tendency at this layer is
dDFe/dt = alpfe * PHI_DUST (source: dust × solubility)
− scav_rat * 86400 * DFe * POC (sink: particle scavenging)
− Q_FE * growth_total (sink: biological uptake)
alpfe controls the surface concentration directly. scav_rat removes iron at a rate proportional to DFe × POC in the same surface layer. With only surface observations of DFe, the optimizer can compensate a too-low alpfe (less source) with a too-low scav_rat (less sink) and reach the same equilibrium DFe — the two parameters are linearly degenerate in the surface budget. v2.6 confirmed this is structural: across GEOTRACES_W ∈ {0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.50, 1.0} and 10 seeds, scav_rat stays at 0.89–0.92 off Carroll regardless of weight.
The physics that breaks the degeneracy is the vertical profile of DFe. alpfe is a surface-only flux (atmospheric dust deposits at the air-sea interface); scav_rat acts throughout the water column wherever POC sinks. The two parameters therefore have different vertical signatures:
alpferaises DFe at the surface but does not directly affect deep DFe.scav_ratsteepens the surface-to-deep DFe gradient (more scavenging → more depletion as POC sinks through the column) and shapes the slope between the ferricline maximum (~500–1000 m) and the deep ocean.
2. Empirical motivation: GEOTRACES depth profile in Eq Pacific AOI¶
From D:\geotraces\GEOTRACES_IDP2025_Seawater.nc, Eq Pacific AOI (lat ∈ [−10, 10], lon ∈ [−180, −130]), Fe_D_CONC with SeaDataNet QC = 49 or 50:
| Depth band (m) | n samples | Fe_D mean ± std (nmol/kg) | Physical interpretation |
|---|---|---|---|
| 0–50 | 67 | 0.162 ± 0.104 | Surface — biological depletion |
| 50–100 | 58 | 0.144 ± 0.095 | Mixed-layer base |
| 100–200 | 93 | 0.196 ± 0.124 | Thermocline gradient |
| 200–500 | 102 | 0.414 ± 0.145 | Regeneration onset |
| 500–1000 | 92 | 0.599 ± 0.115 | Ferricline maximum |
| 1000–2000 | 66 | 0.574 ± 0.103 | Mid-water |
| 2000–6000 | 130 | 0.488 ± 0.089 | Deep (slow ligand-buffered scavenging) |
| Total | 608 |
This is a classic Pacific iron profile (e.g. Tagliabue et al. 2017 Nature; Resing et al. 2015 Nature): low surface DFe (biology), rising through the thermocline as POC sinks and remineralizes, peak at the ferricline (~500–1000 m, where remineralization rate balances scavenging), then slow decline through the deep ocean. The peak-to-deep slope (0.60 → 0.49 nmol/kg over ~3000 m) is the scavenging-dominated regime.
v2.6 uses only the 67 surface samples and bins them onto Darwin's 21×51 Eq Pac grid (25 populated cells). The remaining 541 samples are discarded by depth_max = 50 m in geotraces_loader.bin_to_grid.
3. Minimum-viable design: 2-layer box¶
A 2-layer split is the simplest extension that exposes the scav_rat signature. Going to 3+ layers refines the ferricline resolution but does not add a new degree of freedom over the 2-layer model for the alpfe-vs-scav_rat identifiability question. v2.8 may add layers; v2.7 stays at 2.
3.1 Layer definition¶
- L1 (surface, 0–50 m): photic + mixed layer, matching v2.6's
H_MLD = 50 mexactly so that recoveredalpfeis directly comparable to Carroll's published 0.928 without re-scaling the dust-per-volume convention. Hosts all 5 phytoplankton tracers. Receives atmospheric dust flux at the v2.6 PHI_DUST rate. - L2 (subsurface, 50–1000 m): twilight zone + ferricline. Phytoplankton zeroed (HNLC Eq Pacific has a shallow DCM; subsurface biology is negligible for iron dynamics at scoping precision). Layer thickness
h2 = 950 m. Hosts 297 GEOTRACES samples (50–1000 m bands aggregated).
Below 1000 m is excluded from v2.7. The deep ocean (1000–6000 m, 196 samples) constrains long-residence-time scavenging which is dominated by ligand chemistry, not the scav_rat * DFe * POC term in Carroll's parameterization. Adding a deep box would require introducing ligand physics; v2.7 stays within Carroll-6.
3.2 State vector (15 tracers total)¶
L1: [DFe_1, P_diatom, P_lge, P_syn, P_proLL, P_proHL, POC_1, PIC_1, DIC_1, ALK_1] (10)
L2: [DFe_2, POC_2, PIC_2, DIC_2, ALK_2] ( 5)
Phytoplankton in L2 fixed at 0 (no light, no biology at scoping precision). PIC/DIC/ALK in L2 evolve via remineralization of sinking POC; including them keeps the carbonate system internally consistent even though they don't feed the loss (the carbonate loss is surface-only in v2.7).
3.3 Vertical-coupling physics (NEW vs v2.6)¶
Iron accounting (post-fix v2.7.1). The model does not track Fe-on-POC as a separate tracer, so iron leaves DFe_1 ONLY via biological uptake (fe_uptake = Q_FE × growth_total). The phyto → POC → sinking-POC iron pathway is implicit — iron rides along with POC mass but is not double-bookkept. Iron returns to dissolved phase ONLY in L2, when POC_2 remineralises (fe_remin_L2 = Q_FE × r_remin × POC_2). Whatever fraction of sinking POC does not remineralise in L2 sinks to deep and is lost to the system (consistent with the W_SINK·POC_2 boundary outflow).
The v2.7.0 draft (pre-fix) also subtracted a separate Q_FE × W_SINK × POC_1 term from DFe_1 and credited the same quantity to DFe_2; Greptile correctly flagged this as iron double-counting on PR #42 at e19a5d0. The fix on eccf6d8 → b91f5b8 removes both spurious terms.
Vertical eddy diffusion. Diffusion across the thermocline mixes DFe at rate Kz × (DFe_2 − DFe_1) / dz_int, where dz_int = (h1 + h2) / 2 = 500 m and Kz ≈ 1e-5 m²/s ≈ 0.864 m²/day (typical thermocline value; Whalen et al. 2012). This flux is positive into L1 when DFe_2 > DFe_1 (which the observed Eq Pac profile says is the case). Order-of-magnitude check: with a 0.4 nmol/kg gradient the diffusive flux is ~1000× smaller than the dust source, so it's included for closure but not the dominant coupling.
Remineralisation in L2. Sinking POC remineralises at rate r_remin = 0.05/day (Martin curve b ≈ 0.86; remineralization scale ~200 m in the upper twilight zone). Remineralisation releases the iron implicitly attached to sinking POC back to DFe_2 at the Q_FE = 1e-5 mol Fe / mol C stoichiometry. This is the SOLE iron return path to L2 in v2.7.1.
3.4 Per-layer tendencies (post-fix v2.7.1)¶
L1 iron (surface, 0–50 m):
dDFe_1/dt = alpfe * PHI_DUST [source: dust, per-volume v2.6 convention]
− scav_rat * 86400 * DFe_1 * POC_1 [sink: surface scavenging]
− Q_FE * growth_total [sink: biology uptake — implicitly exports iron to L2 via POC]
+ Kz * (DFe_2 − DFe_1) / (dz_int * h1) [diffusive exchange with L2]
h1 = 50 m matches v2.6 exactly so alpfe recovery is directly comparable to Carroll's 0.928 without re-scaling.
L2 iron (subsurface, 50–1000 m):
dDFe_2/dt = + Q_FE * r_remin * POC_2 [source: POC remineralisation (the implicit iron return)]
− scav_rat * 86400 * DFe_2 * POC_2 [sink: subsurface scavenging]
− Kz * (DFe_2 − DFe_1) / (dz_int * h2) [diffusive exchange with L1]
Conservation: L1 lost via fe_uptake ≈ Q_FE × growth_total; L2 gained via fe_remin_L2 = Q_FE × r_remin × POC_2. At steady state POC_2 = (h1/h2) × W_SINK × POC_1 / (r_remin + W_SINK), so the L2 iron source scales with the iron extracted from L1. The fraction not remineralised in L2 (W_SINK / (r_remin + W_SINK) ≈ 67% at the v2.7 defaults) sinks to deep and is lost.
POC export and remineralisation:
dPOC_1/dt = (biology terms from v2.6) − W_SINK * POC_1
dPOC_2/dt = (h1/h2) * W_SINK * POC_1 − r_remin * POC_2 − W_SINK * POC_2
Same structure for PIC (with R_PICPOC stoichiometry); DIC (with respiration); ALK (carbonate dissolution).
Identifiability claim — what we WERE hoping for. The L2 DFe target was meant to constrain scav_rat independently: with the L2 budget set by inflow (POC remin) balanced against scavenging, GEOTRACES L2 observations would pin scav_rat. In practice (Section 7), the L2 source rate after the fix is much smaller than the pre-fix double-counted version, and scav_rat recovery does not improve enough to break Cal-grade — the structural-limit interpretation from v2.6 still applies.
3.5 Loss function (NEW vs v2.6)¶
Surface terms (unchanged from v2.6): 11-target z-score loss on L1 DFe / Chl1-5 / POC_1 / PIC_1 / DIC_1 / ALK_1 / CO2_flux, PINN drift loss on L1 iron budget, GEOTRACES surface absolute-units MSE on L1 DFe (using 25 populated surface bins from v2.6).
NEW: GEOTRACES subsurface absolute-units MSE on L2 DFe. Aggregates 50–1000 m GEOTRACES samples (~287 samples) onto the Eq Pac 1°×1° grid by bin-mean, then computes
L_geo_subsurface = (1 / N_cells_subsurface) × sum((DFe_2 - GEOTRACES_L2)² × mask) / mean(GEOTRACES_L2²)
normalized by the subsurface mean-square magnitude (same convention as the surface term), with a fresh weight GEOTRACES_W_SUBSURFACE separately tunable. v2.7 default: GEOTRACES_W_SUBSURFACE = GEOTRACES_W = 0.3 (mirror the surface v2.6 winner; sweep TBD if results don't recover scav_rat).
4. Experimental plan¶
- Implementation sanity check. Verify the 2-layer integrator is autograd-clean (gradient flows from
paramsthroughstatetoloss); confirm steady-state DFe profile matches GEOTRACES qualitatively whenparams = Carroll. - Baseline experiment. Run 3 seeds at the v2.6 winner config (
GEOTRACES_W = GEOTRACES_W_SUBSURFACE = 0.3,PINN drift w = 3.0, 1500 epochs, RTX 5090). Compare recovered parameter vector to v2.6 (seed-averaged) and Carroll's published values. - Headline metric. Does
scav_ratmove from 0.92 off → Cal-grade (≤ 0.40)? Doesalpfestay Cal-grade or degrade? DoesR_PICPOCstay drifted (expected — no independent PIC observations were added)? - Subsurface weight sensitivity. If baseline 5/6, declare v2.7 done. If 4/6 (scav_rat still stuck), sweep
GEOTRACES_W_SUBSURFACE ∈ {0.5, 1.0, 2.0}to see if more subsurface anchoring breaks scav_rat free. - Honest reporting. Any weakening of v2.6 recovery (alpfe / Smallgrow / diatomgraz drifting because the loss landscape changed) is reported, not hidden. Multi-layer is only worth the refactor if it net-improves recovery — moving the stuck pair while keeping the recovered triple is the success criterion.
5. What's NOT in v2.7 (deferred)¶
- 3+ layer column model with explicit Martin-curve remineralization and resolved ferricline. The 2-layer model captures the surface-vs-subsurface degeneracy break; finer resolution is v2.8 if needed.
- Depth-resolved POC / PIC / DIC / ALK targets from LLC270 native output. The v2.7 z-score loss stays on surface only; only DFe gets a subsurface target.
- Ligand chemistry for the deep (>1000 m) regime. Outside the Carroll-6 parameter scope.
- Cross-basin validation (N Atlantic, S Ocean). Cluster-gated.
- Independent PIC observations for R_PICPOC unblock. Half-week loader work, separate PR.
6. Expected outcome¶
Success criterion: scav_rat moves to Cal-grade (≤ 0.40 off Carroll) without degrading the v2.6 alpfe/Smallgrow/diatomgraz recovery. Headline becomes "5/6 cal-grade in Eq Pacific" with R_PICPOC as the lone stuck parameter pending independent PIC observations.
Risk: If the L2 DFe target is too sparse on the Eq Pacific grid (the 287 subsurface samples may concentrate at a few stations), the L2 loss may be locally unconstrained and scav_rat may stay stuck. We'll know from the populated-cell count after binning.
Falsification: If scav_rat stays at 0.92 off Carroll even with a clean subsurface DFe constraint, the surface-only degeneracy claim from v2.6 was incomplete — there is another structural source of unidentifiability we haven't named. Worth documenting honestly.
7. Experimental results (11 runs, 5+3+3 seeds across 3 configs)¶
All runs: Eq Pac AOI (lat −5 to 15, lon −160 to −110), 21×51 grid, 1071 ocean cells. Loss config: 11-target z-score (L1 only) + PINN drift weight 3.0 + GEOTRACES surface MSE weight 0.3 (26 populated 1° bins from 52 ≤50 m samples). 1500 epochs, Adam lr 5e-3, ~90 s/run on RTX 5090 Laptop with torch.compile. New v2.7 subsurface loss: GEOTRACES 50–1000 m MSE (28 populated 1° bins from 297 samples) at weight GEOTRACES_SUB_W ∈ {0, 0.3, 1.0}.
7.1 Per-config recovery (post-fix v2.7.1; seed-mean ± std on the recovered value)¶
| Config | n | alpfe | scav_rat (×1e-7) | Smallgrow | Biggrow | diatomgraz | R_PICPOC | Cal-grade |
|---|---|---|---|---|---|---|---|---|
| v2.6 baseline (surface only, 10 seeds, no v2.7) | 10 | 0.549 ± 0.069 (5/10 Cal) | 0.465 ± 0.076 (0/10) | 0.661 ± 0.102 (8/10, 3 Exc) | 0.248 ± 0.056 (3/10) | 0.791 ± 0.137 (9/10) | 0.166 ± 0.010 (0/10) | 3-4/6 |
| SUB_W=0 (2-layer geometry only) | 3 | 0.578 ± 0.031 (2/3 Cal) | 0.451 ± 0.045 (0/3) | 0.729 ± 0.040 (3/3 Cal, 1 Exc) | 0.237 ± 0.040 (1/3 Cal) | 0.746 ± 0.145 (3/3, 2 Exc) | 0.172 ± 0.007 (0/3) | 2-4/6 |
| SUB_W=0.3 (baseline v2.7) | 5 | 0.601 ± 0.027 (4/5 Cal) | 0.572 ± 0.284 (0/5) | 0.646 ± 0.240 (4/5 Cal) | 0.249 ± 0.061 (3/5 Cal) | 0.794 ± 0.144 (5/5, 1 Exc) | 0.147 ± 0.061 (0/5) | 2-4/6 |
| SUB_W=1.0 (strong subsurface anchor) | 3 | 0.767 ± 0.018 (3/3 Cal) | 0.762 ± 0.046 (0/3) | 0.722 ± 0.009 (3/3 Cal) | 0.185 ± 0.027 (0/3 Loose) | 0.615 ± 0.061 (3/3 Cal) | 0.171 ± 0.011 (0/3) | 3/6 |
| Carroll published | — | 0.928 | 6.025 ×1e-7 | 0.661 | 0.431 | 0.830 | 0.0425 | — |
7.2 Headline finding (corrected): v2.7 is comparable to v2.6, not a leap¶
The pre-fix v2.7.0 draft of this doc claimed alpfe reached 3/3 Excellent (mean offset 0.011) at SUB_W=1.0. That headline was an artifact of the iron-double-counting bug Greptile flagged: removing iron from DFe_1 twice (via fe_uptake AND fe_sink_out_L1) inflated the optimizer's pressure to push alpfe high. After the fix, that recovery does not survive. The corrected alpfe numbers are:
| SUB_W | alpfe mean offset (post-fix) | alpfe Cal-grade rate | (pre-fix headline, for reference) |
|---|---|---|---|
| v2.6 (surface only) | 0.408 | 5/10 | — |
| 0.0 (geometry only) | 0.377 | 2/3 | pre-fix said 0.190 (1 Exc) |
| 0.3 (baseline) | 0.353 | 4/5 | pre-fix said 0.121 (1 Exc) |
| 1.0 (strong) | 0.174 | 3/3 | pre-fix said 0.011 (3/3 Exc) |
The 2-layer GEOMETRY still nudges alpfe modestly closer to Carroll's 0.928 than v2.6 (0.41 → 0.18 at strong subsurface weight) but does not push it into Excellent. The Excellent claim is retracted.
7.3 What the fix actually improved: biology parameters¶
The bug's main scientific effect was suppressing Smallgrow and Biggrow recovery. With the fix in place:
| Param | v2.6 (surface) | v2.7.0 pre-fix SUB_W=0.3 | v2.7.1 post-fix SUB_W=0.3 | Direction |
|---|---|---|---|---|
| Smallgrow | 0.001 (8/10 Cal, 3 Exc) | 0.215 (5/5 Cal, 0 Exc) | 0.023 (4/5 Cal mean Excellent) | ✓ restored |
| Biggrow | 0.426 (3/10) | 0.721 (0/5 Loose) | 0.424 (3/5 Cal) | ✓ restored |
| diatomgraz | 0.047 (9/10) | 0.115 (5/5) | 0.044 (5/5, 1 Exc) | ✓ comparable |
Smallgrow's seed-mean offset of 0.023 at SUB_W=0.3 puts it inside the Excellent band on average (though seed std is 0.240, so the rate is 0/5 Excellent per-seed; the mean Excellent reflects a tight clustering near Carroll). Diatomgraz's seed-mean 0.044 is also Excellent. So the post-fix v2.7 effectively reproduces v2.6's biology recovery rather than regressing.
7.4 scav_rat is essentially unchanged by the fix¶
scav_rat drifts LOW across all configs, the same direction v2.6 already exhibits (corrected from the earlier "drifts high" framing in the v2.6 memory note — direct inspection of v2.6 JSONs confirms LOW drift). The fix does not change this:
| SUB_W | scav_rat (×1e-7) post-fix | offset | band |
|---|---|---|---|
| v2.6 | 0.465 | 0.923 | Drifted |
| 0.0 | 0.451 | 0.925 | Drifted |
| 0.3 | 0.572 | 0.905 | Drifted |
| 1.0 | 0.762 | 0.874 | Drifted |
The subsurface loss term still moves scav_rat slightly toward Carroll, but the effect is small. The structural-limit interpretation from v2.6 stands: the L2 budget is dominated by POC_2 dynamics (which is unconstrained in the model), and the optimizer suppresses scav_rat to match observed L2 DFe regardless of weight. v2.8's depth-resolved POC observations remain the realistic unblock.
7.5 Verdict count comparison (post-fix)¶
Aggregate Cal-grade-or-better count out of 6 parameters per seed:
| Source | 4/6 | 3/6 | 2/6 |
|---|---|---|---|
| v2.6 (10 seeds, surface-only) | 4 | 5 | 1 |
| v2.7 SUB_W=0 (3 seeds) | 1 | 1 | 1 |
| v2.7 SUB_W=0.3 (5 seeds) | 2 | 2 | 1 |
| v2.7 SUB_W=1.0 (3 seeds) | 0 | 3 | 0 |
v2.7 (corrected) is roughly comparable to v2.6 on the aggregate verdict count. At SUB_W=0.3 the median is 3-4/6, matching v2.6. SUB_W=1.0 settles uniformly at 3/6 (consistent, but no upside over v2.6). The 2-layer architecture does not deliver a clean win on this metric.
7.6 Honest summary (post-fix)¶
What v2.7 did contribute:
- A vetted 2-layer integrator with corrected iron mass balance (the bug Greptile caught is now fixed and documented).
- Empirical confirmation that the subsurface GEOTRACES loss alone is NOT enough to break
scav_ratCal-grade. The HNLC structural-limit interpretation from v2.6 is reinforced, not falsified. - An infrastructure base for v2.8 (batched seed runner, target cache, 2-layer integrator) that makes the next research lift cheap to iterate.
What v2.7 retracts:
- The pre-fix "3/3 Excellent alpfe at 0.011 off Carroll" headline. That was a bug artifact and is not the corrected result.
- Any claim that the subsurface GEOTRACES observation alone delivers a substantial improvement over v2.6.
What's still open after v2.7:
scav_ratremains structurally drifted (~0.87-0.93 off Carroll across all configs). Cross-basin (non-HNLC) data is still the realistic unblock, gated on cluster compute.R_PICPOCremains drifted (~0.15 vs Carroll 0.042). Needs independent PIC observations.- POC_2 dynamics are unconstrained, which is the most-actionable v2.8 lever.
Recommended v2.8 next steps:
- Depth-resolved POC observations from LLC270 native output as a subsurface POC_2 target. This is the realistic scav_rat unblock.
- Reconsider the 2-layer structure itself: with POC_2 bounded, would scav_rat finally identify? Or is HNLC truly insufficient regardless of model resolution?
- Cross-basin AOI (N Atlantic) for scav_rat constraint — cluster-gated.
- Independent PIC observations for R_PICPOC — half-week loader work, separate PR.
7.7 Data persisted¶
11 result JSONs at scripts/run_v2.7_multilayer_result_seed{0..4}_surf0.3_sub{0,0.3,1.0}_pinn3.0.json, all regenerated after the fix. Aggregator at scripts/analyze_v2.7_results.py. Test coverage at tests/test_carroll6_5pft_2layer.py (9 tests, all passing against the fixed integrator).