Findings — 2026-05-08¶
The DarwinDiff methodology has been demonstrated end-to-end across nine notebooks. This document consolidates the scientific findings — what the evidence supports, what it does not, and where each number comes from.
1. The architectural argument is quantitatively demonstrated¶
The headline scientific claim of the project is that gradient-based parameter learning beats Green's-functions calibration not because it searches better, but because it can express a richer parametric class — per-cell parameter values conditioned on local environment, where Green's-functions produces one global scalar per parameter for the entire ocean.
On a synthetic two-regime setup (Mid-Atlantic-style iron-replete + Pacific-style iron-limited, with biogeochemistry-plausible regional contrasts in alpfe, scav_rat, Smallgrow, Biggrow, diatomgraz, R_PICPOC):
| Approach | Mean per-region rel error | Loss floor |
|---|---|---|
| Global-scalar fit (Green's-functions parametric class) | 20.8 % | 2.6 × 10⁻¹ |
| Per-region MLP conditioning on (SST, dust flux, MLD) | 1.4 % | 1.7 × 10⁻⁴ |
| Improvement factor | 15.2× | three orders of magnitude |
The global-scalar fit's loss plateau is the structural ceiling — the LSE solution averaging the two regimes' truths. No optimisation algorithm restricted to one global scalar set per parameter can improve on it, because the parametric class itself is too narrow for heterogeneous biogeochemistry. The MLP-conditioned approach is in a strictly larger class that subsumes the global-scalar one and recovers each regime's truth essentially perfectly.
2. Carroll's six parameters are recoverable from synthetic observations within single-digit % — when observations are timed correctly¶
Single-regime recovery against the Carroll 2020 published optimum:
| Parameter | Truth (Carroll JAMES) | Recovered | Rel. err |
|---|---|---|---|
R_PICPOC |
4.245 × 10⁻² | 4.224 × 10⁻² | 0.49 % |
Smallgrow |
0.6610 | 0.6618 | 0.13 % |
Biggrow |
0.4315 | 0.4176 | 3.22 % |
diatomgraz |
0.8300 | 0.7941 | 4.32 % |
scav_rat |
6.025 × 10⁻⁷ | 5.644 × 10⁻⁷ | 6.32 % |
alpfe |
0.9283 | 0.8594 | 7.43 % |
End-to-end autograd through 200 forward-Euler steps × 5 coupled tracers (Fe–Ps–Pl–POC–PIC). Single-shot Adam fit on all six knobs simultaneously.
The crucial methodological caveat: an earlier 200-day spin-up reached the iron pool's equilibrium before the snapshots were taken. At equilibrium, alpfe × dust_source ≈ scav_rat × DFe × POC admits a continuous family of (alpfe, scav_rat) pairs that produce the same DFe — a product degeneracy. With a 50-day spin-up that keeps the iron pool in transient, that degeneracy breaks because the rate of approach to equilibrium depends separately on each rate constant. The recovery shifted from "iron pair stuck at sigmoid-midpoint init" to single-digit %. Observation timing matters more than epoch count for identifiability. This is the same principle as the 4-parameter steady-state-vs-trajectory result from earlier work, now demonstrated on the canonical Carroll-6 target.
3. Per-cell parameter prediction works at 2-D scale¶
Extending the scalar fit to a 128 × 128 grid (16 384 cells) where each cell evolves the same five-tracer box model with its own per-cell parameters predicted from environmental covariates:
| Param | Mean rel RMSE across grid |
|---|---|
Smallgrow |
1.91 % |
Biggrow |
4.55 % |
R_PICPOC |
5.23 % |
diatomgraz |
5.75 % |
alpfe |
10.50 % |
scav_rat |
36.13 % |
Five of six parameters recover to under 11 %. The iron scavenging rate sits at 36 %, carrying the same identifiability gap from the 0-D case. Loss converged 1.96 → 0.0317 over 1500 epochs (~ 60× reduction).
The CPU and GPU runs at this scale produce numerically identical recovery (same RNG seed, same Adam updates) and differ only in wall time. Both devices are dispatch-bound at this grid size — every per-step tensor op operates on 16 K-element tensors, so the kernel-launch overhead (~ 5 µs on GPU, ~ 30 µs Python dispatch on CPU) dominates over the actual arithmetic. The GPU only edges CPU by 1.2× here (422 s vs 511 s). This is the threshold scale where the verdict starts shifting toward GPU; meaningful GPU advantage requires moving to grids in the 10⁶-cell range, where memory-bandwidth-bound regimes dominate.
4. The methodology transfers to real ocean observations¶
Loaded GLODAPv2.2016b mapped climatology (NOAA NCEI, public, no authentication), subset to a Mid-Atlantic AOI (30°N–50°N, 60°W–30°W, 582 ocean cells), and ran the per-cell network against the real-data spatial pattern of surface NO₃. The DINN reads SST per cell as its only environmental covariate and produces per-cell Carroll-6 maps; the box model integrates 200 steps from a uniform initial state at every cell; the loss is MSE between the z-scored predicted phytoplankton biomass field and the z-scored negative NO₃ field — both standardised to zero-mean unit-variance over ocean cells, so the loss is a pure spatial-pattern match decoupled from absolute magnitude.
Result: Pearson correlation r = 0.691 (r² = 0.477) between the predicted box-model phytoplankton biomass spatial pattern and the GLODAP −NO₃ field, over 582 ocean cells. The DINN, trained on SST alone, recovered a parameter-mapping that explains ~ 48 % of the variance in the inverse-nitrate spatial structure (r² interpretation; the 0.69 magnitude is the correlation coefficient itself, not variance explained). Loss converged 3.00 → 0.617 over 1500 Adam epochs.
The recovered Carroll-6 maps are smooth functions of SST (because the network is per-cell with one input feature) and lie in physically plausible ranges, no parameter saturating its bound. They do not match Carroll's published optima — and could not, because GLODAP is a real ocean observation product, not output from Carroll's specific simulator. The recovered numbers are best-fit values for this loss against this real-data spatial pattern; the comparison to Carroll's optima requires a future fit against ECCO-Darwin v5 output specifically (where the simulator that produced the data is the simulator the parameters configure).
A first attempt at the loss formulation used 1 / (NO₃ + ε) as the target. GLODAP's mapped NO₃ has small negative values in some oligotrophic cells (mapping artifact at near-zero gridded means), so inversion produced large outliers that corrupted the loss landscape and saturated parameters at bounds (Pearson r = 0.486). Switching to z-scored fields fixed the conditioning. This illustrates a real-data lesson the synthetic notebooks could not surface: gridded BGC products contain mapping artifacts that the loss function has to be robust against.
5. Compute regime: dispatch-bound at threshold scale, memory-bound at production scale¶
A simple per-op cost model — launch_overhead + bytes / bandwidth — captures the regime transition that determines hardware choice:
| Grid | Cells | Per-op tensor | Predicted GPU (5090) speedup vs CPU | Predicted B200 speedup |
|---|---|---|---|---|
| 128² | 16 K | 64 KB | 6× (measured: 1.2× — Python overhead higher than model assumes at this scale) | 6× |
| 1024² | 1 M | 4 MB | 12× | 21× |
| LLC270 (~ 12 M cells, 1/3°) | 12 M | 50 MB | 19× | 90× |
The model overstates GPU advantage at small scales because it treats per-op overhead as a single constant; the actual Python interpreter dispatch cost on tiny tensors is higher and not fully amortised even at 16 K elements. From 1 M elements upward the model becomes accurate as memory bandwidth dominates the per-op cost.
Memory budget for the autograd graph through forward-Euler integration, with optional activation checkpointing (segment of K timesteps with K × n_intermediates × state_bytes per-segment intermediates plus ceil(T/K) × state_bytes boundary states):
| Workload | No checkpoint | Checkpoint @ 5 | Fits 5090 (24 GB) | Fits B200 (192 GB) |
|---|---|---|---|---|
| Notebook-07 scale (128² × 200 steps × 5 tracers) | 0.6 GB | — | yes | yes |
| Mid-Atl AOI patch (1.5 M cells × 30 days × 39 tracers) | 65 GB | 12 GB | yes (with ckpt) | yes |
| LLC270 global × 30 days × daily | 523 GB | 98 GB | no | yes (with ckpt) |
| LLC270 global × 1 yr × daily | 6 364 GB | 214 GB | no | no (single-GPU) |
The single-GPU memory budget closes for global multi-month autograd on B200 with aggressive checkpointing; it does not close for global multi-year autograd, which requires either multi-GPU parallelism or time-window decomposition (compute and accumulate gradients month-by-month rather than building one autograd graph across the full year).
6. Carroll-6 → Darwin 3 namelist mapping is partially verified¶
Carroll 2020's six parameters were tuned in Darwin 1, where they live as named scalars in hardcoded Fortran. The implementation backend for DarwinDiff is Darwin 3 (used by ECCO-Darwin v5), where parameters are exposed via namelists. The mapping from Carroll-JAMES names to Darwin 3 namelist entries, verified against MITgcm-contrib/ecco_darwin/v06/llc270/input_darwin/data.darwin and data.traits:
| Carroll JAMES name | Darwin 3 namelist target | v06/llc270 value | Status |
|---|---|---|---|
alpfe |
alpfe (&DARWIN_PARAMS) |
1.0 | verified |
scav_rat |
scav_tau (&DARWIN_PARAMS) |
1.0 | partial (units / sign convention not yet confirmed) |
Smallgrow |
a_PCmax(1) PicoCyano / a_PCmax(2) PicoEuk |
1.1098 / 1.2104 | verified |
Biggrow |
a_PCmax(5) Diatom |
2.4967 | verified |
diatomgraz |
PALAT(prey=5, predator=?) (data.traits) |
matrix entry, predator index TBD | partial |
R_PICPOC |
a_R_PICPOC array (&DARWIN_TRAIT_PARAMS) |
4.19 × 10⁻² (default), 5.0 × 10⁻² (Cocco) | verified |
Four of six are pinned; the two open items (scav_tau units, PALAT predator index for diatom grazing) require source-code inspection of pkg/darwin/ rather than further notebook work. They should land before any production fit against real ECCO-Darwin v5 output, otherwise the optimiser would target the wrong knobs.
7. What the evidence does not yet support¶
These claims appear plausible but have not been quantitatively demonstrated:
- Recovery of Carroll's published optima from real observations. The synthetic recovery (notebook 05) and the real-data pipeline test (notebook 09) are different experiments. Recovery on real data requires fitting against either ECCO-Darwin v5 output (where Carroll's optima are the simulator's actual parameter values) or vessel observations including dissolved iron.
- A direct head-to-head against Green's-functions on real data. The 15× win in section 1 is on synthetic two-regime data with planted heterogeneity. Repeating that comparison on real GLODAP or ECCO-Darwin v5 fields with a real Green's-functions baseline (one global scalar per parameter, fit by autograd LSE as the closest computational analog) would be the cleanest extension of the synthetic result.
- Iron-pair identifiability on any real observational target. GLODAP has no gridded dissolved-iron field; the iron pair (
alpfe,scav_rat) is structurally unconstrained by the notebook 09 loss. Adding GEOTRACES iron sections (sparse but real) is the next-step deliverable for those two parameters. - Spatial coupling. All notebooks evolve cells independently. Real ocean biogeochemistry has horizontal and vertical advection, mixing, and biological communication via dissolved organic matter. Adding spatial coupling is required for any global fit.
- Multi-year autograd. No notebook has run the autograd graph beyond a 50-day spin-up. The memory budget shows multi-year scale exceeds single-GPU capacity even on B200 with checkpointing; the architectural choice between multi-GPU parallelism and time-window decomposition has been scoped but not implemented.
8. Where this leaves the project¶
The methodology is now established with quantitative metrics at every stage from synthetic recovery to real-data transfer. The next step is the substitution of ECCO-Darwin v5 output for the synthetic and proxy targets — at that point the Carroll-recovery-on-real-data comparison becomes possible, the iron pair becomes constrained by the real Darwin simulator's iron field, and the head-to-head against Green's-functions can be run on actual ECCO-Darwin tracers. Compute requirements for that fit are scoped: AOI patches with checkpointing fit a laptop 5090; global single-month autograd needs B200 with checkpointing; global multi-year autograd needs multi-GPU or time-window decomposition.