# AdaPBS, v8 gdirs, round 7

(Moved verbatim from CLAUDE.md on 2026-11-02. History and reasoning, not instructions.)

## adapbs/: handover for an ADAPBS (melt_f, prcp_fac, temp_bias) recalibration (2026-10-30)

Why the RGI-date state is wrong: equifinality, not a bug.
- The dyn calibration matches the total mass change 2000–2020 (Δmass / RGI
  area, as Hugonnet), which is right.
- With melt_f as the only free parameter, though, a glacier that is too
  small and loses the right kg fits as well as one of the right size.
- OGGM's remedy is the dynamic spinup (an area match at the RGI date). We
  use the plain `dynamic_melt_f_run` without it.
- In `round5/output/sw_pf3.5_mf5_lam12_ys1500`, Δmelt_f anti-correlates with
  A(2003)/A_RGI (Spearman −0.71). Where melt_f rose by ≥ 2 (26% of the
  area), the glaciers end at 0.5–0.65× their RGI area and 0.25–0.34× their
  inversion volume.
- Great Aletsch: melt_f 6.87 -> 14.15 ("full success"); state at 2004
  56.8 km² / 4.0 km³ against 82.2 / 13.2.
- With the static parameters it is 84.9 / 13.4, but dmdtda is −1658
  against −1211 ± 132.

`adapbs/aletsch_forward.py` is a URL-only forward model (1500 -> 2025 on the
v7 `sw_pf3.5_mf5_lam12` gdirs), tested 2026-10-30 in a clean OGGM master
5e9f1332 environment; it reproduces the round5 run exactly. `adapbs/README.md`
holds the instructions for the user's laptop Claude, also online:
https://cluster.klima.uni-bremen.de/~fmaussion/runs/lmr_2025/adapbs/

### AdaPBS on the cluster (2026-10-31)

Handover from the laptop: OGGM branch `fmaussion/oggm@adapbs` (c9b2f7ad,
`oggm/core/adapbs.py`), scripts in `adapbs/laptop/` (`calib_aletsch.py`,
`analyse.py`; edited here, originals `*.laptop_orig`). Cluster wrappers:
- `adapbs/run_calib_slurm.sh <RGI> <tag> [opts]` (one glacier);
- `adapbs/run_batch_slurm.sh <ids> <tag> [opts]` (array);
- `adapbs/summarise.py [out_dir]` (posterior, Neff, z-misfits, flags, cost).

- **Area discretisation** (`area_quantisation.py`): one terminus grid point
  is a median 15% of the area for glaciers < 0.1 km² (2093 glaciers, 78
  km²), ≤ 2.5% above 0.5 km². `--area_k` adds k·dA² to the area variance
  (0.25 used).
- **Speed:**
  - All nodes have the same CPU (EPYC 7542). One run on one core: Aletsch
    5.8 s, Pizzo Tambau 1.8 s; the laptop does Aletsch in 2.4 s.
  - The calibration script adds 40–90% per particle.
  - OGGM hot spot: `FlowlineModel.get_mb` calls `floatyear_to_date` every
    step (only the year is used), and `get_days_of_month` calls it 12× per
    year. A scalar fast path (`adapbs/fastdate.py`, test `bench/bench2.py`)
    gives 0/1.02M differing dates, bit-identical output and 36% faster runs.
    Worth upstreaming.
- **11 test calibrations** (200 particles, MB σ 0.2 × err, area 5%, vol 30%):
  - Aletsch, Rhône, Hintereis and Gorner fit all three targets (|z| < 0.7).
  - Gepatsch's dmdtda is 2.2σ off. Findelen fails (z −2.8 / −3.3).
    Plaine Morte collapses to Neff 1.5. The small glaciers fit poorly.
  - 100 particles give the same Aletsch calibration as 200.
  - **In all 8 large glaciers, prcp_fac goes down (1.6–2.6, from a static
    2.8–4.2) and melt_f goes down:** less mass turnover. It is systematic:
    the transient fit disagrees with the prcp_fac 3.5 preset.
- **Overnight batch** (array 589003): 317 glaciers (all > 2 km² + 100
  small, stratified; `adapbs/overnight_ids.txt`, 65% of the area); 100
  particles, ≤ 6 iterations, area_k 0.25; fastdate on; vanished glaciers
  no longer count as failed. -> `adapbs/out_night/`; then
  `python summarise.py out_night`.
- **Overnight results (2026-10-31):** 317/317 ran, 0 crashes, 600 runs
  each; 180 core-h (median 27 CPU-min per glacier, 2.7 s per run).
  `out_night/summary.csv`.
  - Posterior-median fit, % of area: dmdtda within 1 Hugonnet σ 94%;
    A2003 within 10% 80%; V within 30% 82%; all three 77% (186
    glaciers). Within the run's own tolerances (0.2σ MB): 54% within 1,
    72% within 2.
  - Against round5 ens010 (same glaciers): A2003 within 10% 23% -> 80%,
    V within 30% 41% -> 82%. Aletsch 58 -> 84 km² (obs 82), Gepatsch 6 ->
    17 (17), Rhône 11 -> 16 (16). Findelen is still too small (9.7 vs
    14.3).
  - Misfits are one-sided: too much mass loss *and* too small (95% / 92%
    of the misses), and the volume is low (median V/V_inv 0.83). Small
    glaciers (< 1 km²) fit worst.
  - prcp_fac drops for 96% of the glaciers (median ×0.68), melt_f drops for
    89%, temp_bias rises for 75%. Nothing is at a prior bound.
  - The 1850 area summed over the 317 drops 1770 -> 1639 km²: big glaciers
    grow, small ones no longer overgrow. Check the regional LIA area
    before adopting this.
  - The `neff<50` flag is meaningless at 100 particles (median Neff 30);
    Neff < 10 is 20% of the area.
- **WGMS test, two priors (2026-10-31).** OGGM `fmaussion/oggm@50c97b88`:
  the adapbs code is unchanged, plus upstream #2024 / #2025 speed-ups, so
  no `fastdate.py`. It reproduces the Aletsch reference runs exactly.
  `ADA_EXP` selects the gdirs (`sw_pf3.5_mf5_lam12` default, or
  `sw_pf2.5_mf5_lam12`).
  - Glaciers: `wgms_ids.txt`, 32 glaciers of 11-01 with ≥ 20 WGMS annual
    balances.
  - Runs: arrays 589377 (pf3.5) and 589378 (pf2.5), overnight options, 13
    / 11 CPU-h -> `out_wgms_pf35/`, `out_wgms_pf25/`.
  - `wgms_compare.py` -> `wgms_compare/`: model B = 900·ΔV / mean A
    (hydro years = WGMS years), detrended std and trend on the common years.
  - Fit (dmdtda 1 Hugonnet σ, A 10%, V 30%): pf3.5 88% of the area, pf2.5
    89%.
  - **The posterior follows the prior along a prcp_fac–temp_bias ridge.**
    Both priors drop prcp_fac by ~0.6–0.9 from their own static value
    (median 3.29 -> 2.36; 2.30 -> 1.71). The ~1 K temp_bias difference of
    the static calibrations survives. 1850 area of the 32: 268 vs 264 km².
  - WGMS detrended std, model / obs, median (IQR):

    | run | std ratio | within 20% | trend ratio |
    |---|---|---|---|
    | AdaPBS pf3.5 | 0.99 (0.87–1.10) | 59% | 1.15 |
    | AdaPBS pf2.5 | 0.93 (0.78–1.09) | 53% | 1.15 |
    | round5 pf3.5 | 1.40 | 22% | 1.47 |
    | round5 pf2.5 | 1.17 | 56% | 1.38 |

    - Round5's high melt_f values are too variable (Aletsch melt_f 14:
      2.6×).
    - Within a glacier, among acceptable particles, the std ratio rises
      0.19 per unit prcp_fac (r 0.6). Top particles with prcp_fac ≤ 1.6 are
      too smooth (0.82).
    - So WGMS mildly favours the pf3.5 prior, and it constrains prcp_fac
      only to about ±1 per glacier.
    - Detrended r ≈ 0.6 everywhere (temperature only). The glacier-to-
      glacier pattern of the std is poorly reproduced (r 0.38).

## lmr_gdirs_v8_202610: v7 `sw_pf3.5_mf5_lam12` with LAPrec precip (2026-10-31)

User's plan: static calibration only (AdaPBS later), ens010 temperature,
prcp_fac 3.5 / melt_f 5 (LAPrec is drier; AdaPBS will lower prcp_fac), λ12,
b240, consensus, OGGM 5e9f1332; v7 scripts unchanged. Smoke test OK
(589409–589411). **Full build submitted 2026-10-31: tb / build1 / build2
589415 / 589416 / 589417.** **DONE 2026-10-31**, all COMPLETED.
`lmr_gdirs_v8_202610/compare_v7.py` (vs v7 `sw_pf3.5_mf5_lam12`): same 75
failures (geometry), prcp_fac (3.35 area-weighted, 59% of the area at the
±20% bound), melt_f (median 5.0) and inversion volume (127.7 km³) as v7.
Only temp_bias moves: area-weighted −0.39 vs +0.21 (per glacier −0.60,
P10/P90 −0.85/−0.31). The drier LAPrec goes entirely into temp_bias.
Next (not started): runs from these gdirs; AdaPBS (expected to lower
prcp_fac); NAT precip files with `laprec_generator.nat_precip`.

**AdaPBS sanity check on v8 (2026-11-01).** The 32 WGMS glaciers, same
options as the v7 runs (array 589421, 12 CPU-h). `ADA_TREE=lmr_gdirs_v8_202610
ADA_EXP=pf3.5_mf5_lam12` -> `adapbs/out_wgms_v8_pf35/`; comparison with v7
pf3.5 in `adapbs/wgms_compare_v8/` (`WGMS_OUT` env var of `wgms_compare.py`).
- Forcing check (L3 tars, scratch script): r(WGMS winter balance, gdir
  Oct–Apr P 1901–2023) median 0.65 in v8 vs 0.05 in v7 (25 glaciers).
  Temperature identical to v7. P ratio v8 / v7 median 0.8, the same before
  1850 and over 1900–2023. Saas valley 0.59.
- Fit (all three targets): 93% of the area (v7 88%). None fits only in v7;
  Vernagt and Hohlaub fit only in v8. Neither fits: Allalin, Kesselwand,
  Stubacher Sonnblick, Plattalva, Wurten, Kessjen.
- Posterior medians v8 (v7): prcp_fac 2.60 (2.36), melt_f 4.01 (3.81),
  temp_bias −0.20 (+0.32). Per glacier temp_bias −0.68 vs v7 (the static
  calibration already had −0.63). Nothing at a prior bound. Failed runs
  44 (v7 76).
- 1850 area / volume of the 32: 271 / 31.1 (v7 268 / 30.5).
- WGMS: std ratio 1.14 (IQR 0.98–1.28; v7 0.99), per glacier ×1.11 (LAPrec
  has the observed interannual variability). Trend ratio 1.12 (1.15),
  trend r across glaciers 0.67 (0.63). Detrended r 0.61 (0.59): the annual
  balance is summer dominated.
- HEF temp_bias −2.28 (v7 −0.98; prcp_fac 3.14). Ötztal LAPrec is 0.71× v7.
  Plausible, but large.
- Verdict: no big mistake; v8 is at least as good as v7.

**Regional AdaPBS on v8 (submitted 2026-11-01 ~01:30, user's go).**
- Calibration: array **589480** (80 tasks × 8 cpus, 12 h).
  - Glaciers: `adapbs/region_ids_v8.txt`, the 3817 valid 11-01 glaciers
    (99.7% of the area), sorted by area for round-robin.
  - Same options as the WGMS test, plus `--slim`. Each glacier's netCDF
    keeps every particle and weight, with float32 series from 1800, zlib,
    ~0.5 MB.
  - Output `adapbs/out_region_v8_pf35/`; logs
    `adapbs/_logs/region_v8_pf35_589480_*.out`.
  - `run_batch_slurm.sh` now skips glaciers already done (resubmit to
    resume) and reads the gdirs from disk (cache link +
    /home/data/download).
- Gdirs: L4 on v8 `exps/pf3.5_mf5_lam12`, built by
  `lmr_gdirs_v8_202610/adapbs/` (`build_adapbs.py`, `check_adapbs.py`,
  `submit_adapbs.sh`; README there).
  - The top-weighted particle goes into settings.yml, plus `adapbs_*`
    metadata and `adapbs_particles.nc`. Missing calibrations fall back to
    static (flagged).
  - Smoke test 589562 OK: a rerun from the L4 == the stored particle run,
    exactly.
  - **Full build job 589563, afterok 589480** -> `lmr_gdirs_v8_202610/
    adapbs/pf3.5_mf5_lam12_adapbs/`.
- When done:
  - `sacct -j 589480,589563`;
  - `python adapbs/summarise.py adapbs/out_region_v8_pf35`;
  - count `fallback_static` in `<tree>/summary/*.csv`;
  - regional comparison against round5 / v7: 2003 area / volume, Hugonnet,
    the LIA area (1850, 1800–1875 max).
  - If a calibration task fails or times out, afterok cancels the build.
    Resubmit the calibration (done glaciers are skipped), then
    `./submit_adapbs.sh pf3.5_mf5_lam12 region_v8_pf35 pf3.5_mf5_lam12_adapbs`.
- Next: NAT LAPrec precip files (`laprec_generator.nat_precip`), then the
  NAT / HIST runs from the L4.
- **DONE (checked 2026-11-01 21:20).**
  - Calibration: 80/80 tasks COMPLETED, 3817/3817 glaciers, no failure.
    1146 CPU-h (median 18 CPU-min per glacier, 1.8 s per run); 2.6 GB.
  - L4 build: 40/40 COMPLETED. 3817 `top_particle`, 0 fallback, 75
    no_gdir; 2.5 GB.
  - Parameters (median static -> AdaPBS): melt_f 5.0 -> 4.26, prcp_fac
    3.57 -> 2.80 (area-weighted 2.58), temp_bias −0.36 -> +0.10.
  - Fit, % of area (posterior median): dmdtda within 1 Hugonnet σ 95%;
    A2003 within 10% 79% (round5 ens010 20%); V within 30% 77%; all three
    73%. By size, all three: < 0.1 km² 13% of the glaciers, 0.1–0.5 43%,
    1–5 69%, 5–20 85%, > 20 100%.
  - Regional sums: A2003 2014 vs RGI 2083 (round5 2245); V2003 110 vs
    inv 128; dmdtda −843 vs −805.
  - **But the LIA state shrinks** (top particle, the set in the gdirs;
    `adapbs/out_region_v8_pf35_topparticle_lia.csv`):
    - 1800–1875 max area 3389 km² (round5 ens010 3928; round5 HIST
      members 3391–4428). With the +410 km² not-modelled outlines that is
      ~3800, against Reinthaler 4248 ± 217, i.e. ~−11% (round5 +1.8%).
    - **Superseded, like-for-like (2026-11-01, `round7/lia_max_ens010_adapbs.py`
      -> `round7/figs/lia_max_ens010_adapbs.{png,csv}`):**
      - The particle files store OGGM's `area` (any ice thickness), which
        has single-year spikes from thin snow patches and inflates the
        maxima. Round5 / validate.py use `area_min_h` (> 2 m).
      - With round7's HIST ens010 (== the particle run, with area_min_h)
        and validate.py's additive correction (− remnant glaciers + 410
        km² / 9.9 km³ of not-modelled outlines):
        - area 3698 km² (round5 4323; obs 4248 ± 217), i.e. **−13%**;
        - volume 200 km³ (round5 244; obs 285 ± 44).
      - The AdaPBS area target also used `area` (unthresholded): to keep
        in mind for small glaciers.
    - Volume at the max 190 km³ (round5 235); dV LIA -> 2015 −102 km³
      (Reinthaler −180 ± 39; round5 −137). Area ratio LIA / 2015: 1.98.
    - The drop is all in glaciers < 5 km² (LIA area 2338 vs 2845). Glaciers
      > 20 km² are unchanged (243). Glaciers < 0.1 km² are still 6.7×
      their RGI area at the LIA (round5 8.5×).
  - User (2026-11-01): "LIA underestimation is what it is"; maybe better
    with the snow model later, not a priority. Later, the not-modelled
    *volume* (glaciers that disappeared) must also be corrected, not only
    the area.

## round7: ensemble from the AdaPBS gdirs (2026-11-01)

See [round7/README.md](round7/README.md).
- LAPrec forcing for all members:
  - `data/build_laprec_precip.py` (HIST, 192 members) and the new
    `data/build_laprec_nat.py` (NAT, 197 members, into `<laprec>/nat/` so
    that `process_mira_data` reads them unchanged).
  - Jobs 589750 / 589751; ~104 GB.
  - Checked on 17 members: NAT == HIST < 1850, and HIST == ens010 from
    1850.
- Ensemble `ens_v8_adapbs`: 192 HIST + 197 NAT from the L4, runs from
  **1500** (round5: 1000).
  - Smoke test 589752: HIST ens010 == the stored AdaPBS particle run (1e-7),
    and NAT == HIST < 1850.
  - **Full run array 589753** (800 tasks).
- When done:
  - merge on the login node: `MERGE_JOBS=8 MAIN_MEMBER=010 $LOCAL_PY
    round7/merge_ens.py round7/output_ens/ens_v8_adapbs 40`;
  - then `round7/preprocess_ens.py ens_v8_adapbs`;
  - then redo the analysis. The round6 conditional set-up needs new 1850
    states (L5 on the v8 L4), the permutation runs, the HIST L5 runs and
    the commitment runs, adapted from `lmr_gdirs_v7_202610/states/` and
    `round6/`.

### Second AdaPBS round: area 5%, volume 10%, area_min_h (queued 2026-11-01 ~23:00)

User's go after a reweighting test of the stored particles (scratchpad
`reweight/`).
- Reweighting results, sums over the region:
  - Loosening the area tolerance makes things worse: A 10% / V 30% gives
    area 1960, V 108.
  - Tightening the volume helps: A 5% / V 10% gives A2003 2086 (RGI 2083),
    V 117 (inv 128, was 110), area within 10% for 84% of the area,
    Hugonnet within 1σ 91% (was 95%), LIA area +100 km² (~+3%).
  - The Neff of the reweighted particles is low (~12), so a real run
    should do better.
- New option `calib_aletsch.py --area_var min_h`: the area target and the
  stored series use `area_min_h` (> 2 m), not OGGM's `area`. Backup of the
  script: `calib_aletsch.py.bak_20261101`.
- **The chain** (afterok all the way; the round7 ensemble 589753 on the
  first calibration keeps running, merged by 590381):

  | step | job | output |
  |---|---|---|
  | smoke test (1 glacier, 20 particles) | 590379 | `adapbs/out_test_minh/` |
  | gate `adapbs/check_test_minh.py` | 590456 | log `adapbs/_logs/check_minh_*.out` |
  | regional calibration (tag `region_v8_a5v10`) | 590457 | `adapbs/out_region_v8_a5v10/` |
  | L4 build | 590458 | `lmr_gdirs_v8_202610/adapbs/pf3.5_mf5_lam12_adapbs_a5v10/` |
  | round7 ensemble | 590459 | `round7/output_ens/ens_v8_adapbs_a5v10/` |
  | merge + preprocess (`round7/merge_slurm.sh`) | 590460 | `round7/filtered_ens/ens_v8_adapbs_a5v10/` |

- In the morning:
  - `sacct -j 589753,590381,590379,590456,590457,590458,590459,590460`;
  - if a link failed, the rest is cancelled (kill-on-invalid-dep), and
    the calibration resumes on resubmission;
  - then redo `round7/lia_max_ens010_adapbs.py` with the new tree
    (ens010 from `ens_v8_adapbs_a5v10`) and the regional comparisons.
- **DONE (all jobs COMPLETED; chain finished 2026-11-02 03:07):**
  - Both ensembles have 389/389 members merged (9.6 GB each) and their
    `round7/filtered_ens/<name>/` CSVs.
  - Round 2: 3817/3817 calibrated, 0 fallback.
  - `round7/lia_max_ens010_adapbs.py` now reads the merged files and
    compares round5 / AdaPBS 1 / AdaPBS 2. On 3788 common glaciers:

    | ens010 | round5 | AdaPBS 1 | AdaPBS 2 | obs |
    |---|---|---|---|---|
    | area RGI date | 2117 | 2004 | **2083** | 2082 |
    | vol RGI date | 117 | 110 | **118** | 128 |
    | dV 2000–2020 | −47 | −39 | **−40** | −37 ± 5 |
    | LIA area corr. | 4323 | 3698 | **3797** | 4248 ± 217 |
    | LIA vol corr. | 244 | 200 | **208** | 285 ± 44 |
    | dmdtda 2000–2020 (m w.e./yr, per RGI area) | −1.02 | −0.84 | **−0.85** | −0.81 ± 0.10 (Hugonnet) |
    | dA 2003–2015 | −20.6% | −15.6% | **−14.2%** | −14.0 ± 2.9% (Paul et al. 2020 ESSD: 2100 -> 1806 ± 60 km²) |

    The figure has three columns: time series, present day (dA vs Paul,
    MB vs Hugonnet), and the LIA vs the outlines. No insets anymore
    (user, 2026-11-02).
- **`round7/plot_lia_max.py`** (NAME, default `ens_v8_adapbs_a5v10`;
  2026-11-02): the round6 talk figure for round 7 ->
  `round7/figs/lia_max_{area,volume}_<NAME>.png` and
  `lia_max_corrected_<NAME>.csv`.
  - The HIST members are the new 1850 NAT volume quantiles, plus ens010:
    min ens023 (166 km³), P25 ens007, median ens083, ens010 (P59), P75
    ens169, max ens118 (221). These are not round6's members: ens178 is
    now P65, ens123 P1.
  - The right panel is outline-corrected for every member (validate.py
    additive correction, with the remnant ice per member from the merged
    files, cached in `figs/lia_max_remnants_<NAME>.csv`).
  - Results:
    - Area: NAT median 3824 [P5–P95 3616–4121], HIST 3682–4292 (ens010
      3797). 27/197 NAT and 24/192 HIST are within the outlines' 1σ (4248
      ± 217).
    - Volume: NAT 210 [199–225], HIST 200–236; none within 285 ± 44.
  - **Plots start in 1650** (`X0`, user 2026-11-02): the runs start in
    1500 from the RGI geometry, and 1500–1600 is spin-up.
  - **5 highlighted HIST members** (user, 2026-11-02: "5 is a good
    number"): the 1850 quantiles ens023 / 007 / 083 / 169 / 118. ens010
    is only added with `INCLUDE_010=1`. The key member is the top one
    (ens118, P99: corrected LIA area 4292, on the outlines). All 192 HIST
    and 197 NAT members are already run (round 7, both AdaPBS rounds).
- **`round7/attribution/attr_fraction.py`** (NAME, REF_YEAR 2025; 2026-11-02):
  the round 7 version of round6's `attr_fraction_conditional_..._ref2025`
  -> `round7/attribution/figs/attr_fraction_<NAME>_ref2025.{png,csv}`.
  - **There are no permutation runs for round 7 yet.** The NAT side is the
    full round 7 NAT ensemble, the same for the 5 states (the 1850
    quantile members).
  - This was validated on round5 / round6: A against the unconditional
    NAT == A against the conditional one to < 1.5 km³ from 2000 on (< 0.5
    in 2015 / 2025), but up to 6 km³ apart in 1900. So it is valid from
    ~1950 on.
  - Missing outlines: Reinthaler's ΔV minus the not-modelled outlines'
    LIA volume, 180 − 9.9 = 170.1 ± 39, anchored on the model HIST 2015.
  - Results, ens_v8_adapbs_a5v10, 2025:
    - A = 103 km³ (round6 126; NAT median 174 vs 196, HIST 70.6);
    - share per state 108 / 94 / 88 / 83 / 68% (min -> max 1850 state);
    - **Reinthaler + NAT: 53% (16–84% 41–69, 5–95% 35–84)**, against 50%
      without the outline correction and round6's 60% (46–79).
- **The attribution is always conditional from now on (user, 2026-11-02):**
  each 1850 state is restarted with every NAT climate (permutation runs),
  never the full NAT ensemble. With 5 states that is 985 runs; all 192 would
  be ~38,000, which is why the 5 quantile members exist.
  - State tree: `lmr_gdirs_v8_202610/states/a5v10_1850/`, an L5 on the
    AdaPBS round 2 L4. Per state: HIST climate, 1850 glaciers (from 1500),
    and the HIST run from them with geometry.
  - Permutation runs: `round7/perm/`; see the round7 README.
  - Smoke tests exact (state / HIST / twin runs == the ensemble, 0.0).
  - **Chain:** states 591621 -> perm 591622 (400 tasks) -> post 591623
    (merge + collect_hist + preprocess_perm). Then rerun
    `round7/attribution/attr_fraction.py`: it uses `perm_a5v10` when
    present; the unconditional preview is `_uncond`.
- **LIA figures (`round7/plot_lia_max.py`, user 2026-11-02):** the NAT runs
  are no longer shown (after 1850 they are a climate no comparison uses,
  and misleading). The faint lines and the box are all 192 HIST members
  (== NAT before 1850). All-HIST LIA area median 3813 [P5–P95 3599–4124];
  volume 209 [197–225].
- **Shifted-NAT sensitivity (user, 2026-11-02).** Hypothesis: small shifts of the NAT
  baseline after 1850 move the %ANTH more than anything else. `data/build_nat_shift.py <d>`
  -> `data/unstacked/lmr_seas_fullsplit_oldBE_newnat_shift-<d>/` (NAT T − d × linear ramp
  1850->1900, bit-identical before 1850; d = 0.2 and 0.5 K cooler). Figures
  `climate_analysis/nat_spread_alps.py` -> `figs/nat_spread_alps.png`, `nat_shift_alps*.png`
  (0.2 K = 1.0×, 0.5 K = 2.5× the post-1950 member std of the 31-yr mean, 0.20 K).
  Permutation runs on the same states (`NAT_T_URL=<shifted> ./submit_perm.sh a5v10_1850
  perm_a5v10_shift-<d>`; `NAT_T_URL` is now overridable in round7/common.sh): shift-0.2 run
  591970 / post 591971, shift-0.5 run 591972 / post 591973 (queued behind 591622). Then
  `PERM=perm_a5v10_shift-<d> python round7/attribution/attr_fraction.py` and compare with
  perm_a5v10. HIST is unchanged (same L5 HIST runs).
- **For the paper (user, 2026-11-02):** the final, expensive simulations
  will start earlier (e.g. year 1000). The 1000–2000 period is wanted as
  context for the LIA. The 1500 start of round 7 is for development only.


### Conditional attribution, perm_a5v10 (DONE 2026-11-02 ~16:45)

- Runs 591622 (400/400 OK). The slurm post job 591623 sat in the queue (priority), so it
  was cancelled and run on the login node (merge 985 runs + collect_hist +
  preprocess_perm, ~20 min, no NaN). Log `round7/perm/_logs/perm_a5v10/post_login.out`.
- `attr_fraction.py` -> `round7/attribution/figs/attr_fraction_ens_v8_adapbs_a5v10_ref2025.{png,csv}`:
  2025 A = 103.5 km³ (states 101.5–107.4; NAT 16–84% 85.9–122.3), HIST 70.6 km³.
  Share per state (ens023 / 007 / 083 / 169 / 118): 106 / 92 / 87 / 85 / 72%.
  **Reinthaler + NAT: 52% (16–84% 41–69, 5–95% 34–84)**; the unconditional preview gave
  53%. As in round 6, A hardly depends on the state (the state's relaxation is in both
  HIST and NAT).
- Shifted NAT: shift-0.2 (591970) running, shift-0.5 (591972) queued. Merge their outputs
  on the login node too (the post jobs 591971 / 591973 will wait on priority).
- **Key figures redone for round 7 (2026-11-02):** round6's `regional_conditional.py`,
  `conditional_temperature.py`, `explainer_conditional.py`, `explainer_two_states.py`
  copied to `round7/attribution/` (perm_a5v10, states 023/007/083/169/118; default state
  ens083 = median, pair 023 / 118 = extremes; temperature from
  `climate_analysis/series/nat_jjas_glacw_newnat_r7.csv`). Outputs in
  `round7/attribution/figs/`: `attr_temperature_vs_glacier_conditional_perm_a5v10_ens083_w15[_length[_smb]].png`,
  `explainer_conditional_perm_a5v10_ens083_w15.png`,
  `explainer_two_states_perm_a5v10_ens023_ens118.png`, `regional_*_conditional_*`.
  Conditional SNR (round 6 in brackets): volume ToE 1998 (1997), SNR 2021 4.46 (4.74);
  area ToE 2000 (2000), 3.46 (4.30); length 2000 (1999), 4.09 (4.91); SMB ToE 1994 (1994),
  5.21 (4.24); temperature ToE 1988, 8.04 (same forcing). ToE is again the same for all 5
  states (1998 volume, 2000 area); the old method gave 1957–2009. Gap natural − historical
  2015: ens023 74, ens118 79 km³. Unlike round 6, the natural medians of the two extreme
  states do not fully converge (~170 vs ~178 km³ in 2021).
- **Why the natural medians of the extreme states don't converge (2026-11-02,
  `round7/attribution/nat_convergence.py` -> `figs/nat_convergence_glaciers.png`,
  `figs/nat_convergence_per_glacier_state118_vs_023.csv`):** no missing glaciers (0 of 990
  runs). Paired (same NAT climate) state 118 − state 023 = 54.9 km³ in 1850, 18.8 in 1900,
  10.3 in 1950, 8.1 in 2005 (083 − 023: 0.4 in 2005; HIST 118 − 023: 3.8). Of the 8.1:
  - 161 small glaciers (21.7 km², median 0.04 km², ~300 m vertical extent, mean elevation
    2840 m vs 3030 m for the rest): **no decay, two stable sizes under the same climate**.
    Their difference grows (3.4 -> 5.1 km³) in all 197 climates: from 118 they keep growing
    (median 1.34× their 1850 volume by 2005), from 023 they stay small or vanish (58 at
    < 10% of the 118 volume). Height–mass-balance feedback on low-relief glaciers near the
    ELA. Largest: RGI60-11.02249 (2.9 km²), 00632, 00752, 02105, 02054.
  - 142 slow-decay glaciers (194 km², e.g. RGI60-11.01328, 24 km²: 0.94 -> 0.11): 2.0 km³.
  - the other 3511: 1.0 km³.
  Not checked: whether round 6 had the same (its medians converged).
- **Share with 2015 as reference year (2026-11-02, `REF_YEAR=2015 python attr_fraction.py`
  -> `figs/attr_fraction_ens_v8_adapbs_a5v10_ref2015.*`):** as round 6's main figure, only
  the NAT runs that reach 2015 (no held values). A = 76.4 km³ (states 74.5–79.8; NAT 16–84%
  58–95), HIST 98.7 km³; per state 110 / 91 / 84 / 82 / 66%; **Reinthaler + NAT 45%
  (16–84% 32–63, 5–95% 26–80)**, 42% without the outline correction (round 6: 55%,
  40–76). The 2025 version (52%) needs every NAT run held at its last-5-yr mean after it
  ends (2005–2020). Both kept.
- **Shift-0.2 NAT (DONE 2026-11-02 ~17:45, login-node post via `round7/perm/post_login.sh`):**
  2025 A = 128.1 km³ (states 125.2–133.1) vs 103.5 unshifted, i.e. +24.6 km³ for a 0.2 K
  cooler NAT; **Reinthaler + NAT share 65% (16–84% 51–84)** vs 52%. 4 of 990 runs have one
  glacier with NaN. Shift-0.5 (591972) post-processing queued the same way.
- **Per-glacier, round 7 (2026-11-02, `round7/glaciers/`, README there):** ens118, 2025.
  Detection (SNR > 2 for good): volume median 2001 (area-weighted 1997), 81% of the
  glaciers / 99% of the area by 2021; area 2006; length 2006 (ens083 the same years).
  Share 1850–2025: per-glacier median 58% (16–84% 41–74), loss-weighted = regional 72%;
  Reinthaler-emulated (k = 1.42) median 42% (29–55), regional 54%; Spearman with size
  0.70. Y100 median 1947 (area-weighted 1942, regional 1953; round 6: 1934 / 1943).
  A per glacier: ens118 vs ens083 r = 0.998 (regional 107.3 vs 101.2 km³).
- **Validation, round 7:**
  - Leclercq relative length (`climate_analysis/leclercq_alps/compare_oggm_length_1850.py`,
    new `TAG` env var; Bossons fixed to RGI60-11.03646): median r (normalised, at obs
    years) 0.88 ens010 / 0.90 ens118 (round 5 ens010: 0.87); range model / obs 0.57 / 0.65
    (round 5 0.73): the model retreats too little, more so than round 5.
    -> `figs/oggm_vs_obs_length_1850_norm_r7_ens{010,118}.png`.
  - WGMS (`WGMS_OUT=wgms_compare_v8_a5v10 python adapbs/wgms_compare.py region_v8_a5v10`;
    top particle == round 7 HIST ens010): detrended std ratio 1.16 (IQR 0.99–1.37; 50%
    within 20%), trend ratio 1.09, detrended r 0.61 (round 5: 1.40 / 1.47 / 0.59).
- **Shift comparison figure (2026-11-02):** `round7/attribution/nat_shift_effect.py` ->
  `figs/nat_shift_effect.{png,csv}`: state ens083, HIST vs natural median (16–84%) as
  delivered and 0.2 K cooler (0.5 K added when available), A(t) to 2021, share 2025 vs the
  shift. A 2021 (ens083) 88 -> 111 km³; share 2025 52 -> 65%.
- **Photo pairs, round 7 (2026-11-02, `round7/glaciers/photo_pairs.py`, ens118):** Gauli
  (RGI60-11.01275, 2 undated pairs) is left out (missing from one run, so not in the common
  subset). 199 dated pairs, 77 glaciers. Area share between the photos: median 90% (16–84%
  72–103; round 6 94%); 185 of 199 beyond natural variability (round 6: 191). By "before"
  period: ≤ 1900 77%, 1901–1930 84%, 1931–1960 94%, 1961–1990 106%, after 1990 101% (round
  6: 82 / 90 / 99 / 108 / 100). Lowest again the late-relaxing glaciers: Gosau, Hochkönig,
  Hallstätter, Unteraar, Pasterze, Albigna, Gries. -> `round7/glaciers/figs/photo_pairs_*`.
- **Y100 figure, round 7:** `round7/glaciers/figs/attr_y100_perm_a5v10_state118_ref2025.png`
  (made by attribution.py). Regional 1953 (round 6 ens010 1940, ens178 1943), median 1947.
  More glaciers at the ends than round 6: Y100 = 1850 (natural never below its 2025 level)
  127 glaciers / 0.7% of the area (round 6 ens178: 4), 121 of them the bistable "no decay"
  glaciers of nat_convergence.py (they grow naturally from the large ens118 state);
  Y100 = 2025 (natural still losing at the end) 629 / 8.5% (round 6: 243 / 5.7%), e.g.
  Unteraar 11.01328, 11.02801, 11.00719.
- **Shift-0.5 result (2026-11-02, perm 591972, post on the login node):** 990 runs, subset
  3814 glaciers. A 2025 = 169.9 km³ (states 166.5–174.7), HIST 70.6. Share 2025 (Reinthaler
  + NAT): 86% (16–84% 68–111). The natural runs gain ice: median 237 km³ in 2021 for
  ens083 (1850: 187), so per-state shares exceed 100% for every state except ens118.
  Summary 0 / 0.2 / 0.5 K cooler: A 2021 (ens083) 88 / 111 / 154 km³, share 52 / 65 / 86%,
  i.e. ~+7 points per 0.1 K, roughly linear. 0.2 K is 1× the across-member std of the
  31-yr NAT mean after 1950. -> `round7/attribution/figs/nat_shift_effect.{png,csv}`.
- **SNR smoothing test (2026-11-02, `round7/attribution/snr_smoothing.py`, perm_a5v10):** the
  code's 15-yr bar is a trailing window (pandas rolling), the record ends 2021, and the NAT
  members drop out as their climate ends (197 in 2005 -> 127 from 2015). ToE at SNR > 2,
  the same for all 5 states (±1 yr): volume (a) trailing 15 on N and H 1998, (b) raw 1992,
  (c) N smoothed only 1992; area 2000 / 1993 / 1993. SNR 2020 volume 4.3 / 4.9 / 5.2.
  H's interannual noise is 0.14 (volume) and 0.09 (area) of the NAT spread. Temperature:
  (d) trailing 15 1988, (e) centred 31 1977 (SNR ends 2006), (f) trailing 31 1992; SNR 2020
  7.6 / – / 8.7. Glacier minus temperature ToE: now (a−d) +10 (volume) / +12 (area);
  proposed (b−f) 0 / +1. The gap depends on the threshold: SNR > 1 −37, 1.65 0,
  2.58 +3 / +10. The glacier SNR sits on a plateau around 1.7 from 1945 to 1990.
  -> `figs/snr_smoothing_perm_a5v10{.png,.csv,_timeseries.csv}`. Doc not changed yet.
- **Key figure, new definitions (2026-11-02, `round7/attribution/conditional_temperature_newdef.py`):**
  glacier annual, temperature and SMB 31-yr trailing, 15-yr dashed. ToE (SNR > 2, ens083):
  temperature 1992 (15-yr 1988), volume 1992 (1998), area 1993 (2000), length 1994 (2000),
  SMB 2004 (1994). -> `figs/attr_temperature_vs_glacier_conditional_perm_a5v10_ens083_newdef_length_smb.{png,csv}`.
- **Question A vs B, model only (2026-11-03, `round7/attribution/question_a_vs_b.py`, doc
  comms/ v261103):**
  - A: every HIST change vs every NAT change, each from its own 1820–1850 mean
    (HIST i == NAT i before 1850). 2025 (NAT held): NAT median ΔV −10 km³, HIST −114,
    f_A = 92% (pairs 16–84% 76–109). 2015: 89% (69–112).
  - B: model-only shares (X_S − H as the denominator, no Reinthaler): 2025 ens023 106 /
    007 92 / 083 87 / 169 85 / 118 72%. 2015: 109 / 90 /
    83 / 81 / 66. The anthropogenic loss is the same in A and B (~104 km³);
    only the denominator differs.
  - The natural median is not flat: −10 km³ from 1820–1850 to 2025, so f_A < 100%.
  - The 52% headline uses Reinthaler as the denominator.
  -> `figs/question_a_vs_b_ref{2025,2015}.{png,csv}`.