# Rounds 3–4, v5 gdirs (Moved verbatim from CLAUDE.md on 2026-11-02. History and reasoning, not instructions.) ## Round 3: `lmr_seas_fullsplit_oldBE` HIST realisations Runs live in `round3/` (forcing read via `round3/mira_forcing.py`, which downloads the unstacked files from the cluster web server mirroring `data/unstacked/`). Outputs in `round3/output/lmrseas_oldBE/`, aggregated CSVs in `round3/filtered/lmrseas_oldBE/`. - `Tp_monthly.nc` (+ `P_monthly_BE.nc`) — the original single HIST: the previously chosen random member, blended with obs. Its 1850 glacier state was low (~6th percentile of the NAT ensemble), hence the alternative realisations below. **Deleted on 2026-10-04** (orig and unstacked `tp_monthly.nc` / `p_monthly.nc`), at the user's request. Its pre-1850 temperature matched none of the 197 NAT members. - Consequence: the v3/v4 gdirs and round3's main hist run can no longer be rebuilt; their outputs remain. - The v4/v5 `lmr_forcing.py` path without an `ens_member` no longer works. - `TpNAT_monthly_ens###.nc` — 197 NAT LMR members (numbers have gaps). - Alternative HIST realisations = NAT member ### before 1850, blended with BE from 1850 on. Members chosen by their NAT-run 1850 reg-11 volume (from `volume_km3_nat_ens_reg11.csv`): | 1850 vol. quantile | km³ | member | |---|---|---| | min | 212.90 | ens123 | | 0.25 | 235.14 | ens053 | | median | 242.35 | ens183 | | 0.75 | 249.61 | ens153 | | max | 273.32 | ens178 | ens124 is also delivered but is not in this set (a random extra, ~10th percentile, 228.99 km³). - `monthly_hists_v20261001/` — **current** version (delivered 2026-10-01 as `output_LMRens###.mat`, converted to `Tp_monthly_hist_ens###.nc`, then unstacked). Verified: identical to the matching `TpNAT_monthly_ens###` for all years < 1850. **The blend is gradual over 1850–1899, not a switch at 1850** (checked 2026-10-04 at 10.5°E 46.5°N). The HIST members are identical to each other only from 1900. In 1850–1899 they still differ from each other, by a decadal-mean |Δ| of ~1.05 K in the 1850s -> 0.13 K in the 1890s -> 0 in the 1900s, and also from their NAT member. So the HIST glacier runs can diverge strongly in 1850–1950 (in lam8_dyn, ens153/183 grow to ~250 km³ around 1870–1890 while ens178 drops to ~165 km³). That comes from the forcing, not a model problem. Checked over the full grid, all months: the max |member − ens183| temperature decreases linearly, 2.41 / 1.93 / 1.45 / 0.96 / 0.48 K in the 1850s–1890s, i.e. a 50-yr linear ramp from the member's LMR to BE; it is exactly 0 from 1900. **Precipitation is never shared:** every member (HIST and NAT) uses its own `ensemble/p_monthly_ens###.nc` over the whole period, with no obs blending. Decadal means differ by up to 17–27% at some cells, also after 1900; glacier-area-weighted they differ by only ~±3%. **The 1850–1900 transition is by design:** the provider's own `data/orig/lmr_seas_fullsplit_oldBE/validation_timeseries_JJAS_alpswide.png` marks 1850 and 1900 with dashed lines. Berkeley Earth itself goes back to 1750 (land), so "BE starts in 1900" is not the reason; most likely the obs are only fully trusted from 1900 (inference, unconfirmed). The member-specific part in 1850–1899 is *not* the member's NAT series: HIST = NAT up to and including 1850, but the HIST member differences in 1855–1899 are not proportional to their NAT differences. Probably an all-forcing LMR series that was not delivered. Berkeley Earth has no precipitation product, so all precip is LMR-based. (`P_monthly_BE.nc`, now deleted, was the old random member's precip, not BE obs.) Open question for the provider: what is the 1850–1899 member part? The 1850–1950 divergence between HIST members is **mainly temperature**. Glacier-area-weighted, in 1850–1874 ens153/178 are 5.33/5.12 °C against 3.98/3.99 for ens053/183 (~1.2 K warmer), and their volume drops to 164/180 km³ by 1880 while ens053/183 grow to 236/247 km³. The precip differences are small and do not line up with the volume differences. - **`monthly_hists_v20261005/` and `jjas_hists_v20261005/`** — Mira's fix (delivered 2026-10-05 as `Tp_monthly_LMRens_###[_monthlyBE].nc`, 6 members incl. ens123; `unstack.py`'s `ENS_RE` now also matches this name pattern). The bug was a climatological offset when merging HIST + BE. Checked 2026-10-05: - Both new versions are identical to the NAT member before 1850, and the members are identical to each other from 1900 (as before). - The 1850–1874 member offset is fixed: glacier-weighted means are 4.62–4.69 °C (monthly) / 4.51–4.71 (JJAS), against 3.98–5.33 in v20261001. - **NEW JJAS == v20261001 exactly from 1900 on** (max |diff| 0.0000 K); only 1850–1899 changed. So v20261001 was the JJAS variant, and the current gdirs and static calibrations (2000–2020) stay valid for JJAS. - **NEW monthly differs after 1900:** 1900–1999 ≈ +0.01 K annual, but 2000–2024 −0.24 K annual (JJA +0.21, DJF −0.49, glacier-weighted). ~7% of grid-months differ by > 3 K (max 12 K): monthly BE anomalies give noisier, colder winters and warmer summers, as Mira warned. Using it needs everything rebuilt from the gdirs (temp bias run, static and dynamic calibration). - `alternative_hists/` — **old, wrong** version, kept for reference only: its pre-1850 data does *not* match the NAT member of the same number (its ens124 is just a copy of `Tp_monthly.nc`). The round3 OGGM `lmr_mira_hist_ens###` outputs produced from it are therefore invalid. `round3/run_lmr_mira_seas_oldEB_hist_ens.py` still points `temp_server` at `alternative_hists/` until switched to `monthly_hists_v20261001/`. ## lmr_gdirs_v5_202610: sensitivity gdirs (ens183 HIST) Region 11 L3 gdirs (border 160, L4/L5 not needed) forced by `lmr_seas_fullsplit_oldBE` HIST ens183 (`monthly_hists_v20261001`), as deltas from a default (prcp_fac 2.5, melt_f 5): prcp_fac 1.5/3.5, melt_f 3/7/9, centerlines, trapezoid_lambdas 6 (from L1). One folder per experiment in `exps/`. **Every experiment has its own temp bias run**, made from the same params.cfg as its build: never share a temp bias file between configs. Chunked prepro on OGGM master + `oggm_informed_fixed_prcp_fac.patch`. See [lmr_gdirs_v5_202610/README.md](lmr_gdirs_v5_202610/README.md). The first 8 experiments are below. Five more were added on 2026-10-03, all elev bands from L1: `lam8`, `lam6_ib`, `lam8_ib`, `lam8_mf7`, `lam8_ib_mf7` (see the round4 section). `ib` = `INV_VOLUME=iceboost` in config.sh; the default is consensus, set in `common.sh`. ### Status (2026-10-03, 16:30) - **Done.** All 8 experiments completed (jobs 573888–573911, 48/48 OK): L3 gdirs in `exps/*/RGI62/b_160/L3`. Area-weighted over 11-01: | exp | failed | temp bias file mean | prcp_fac | melt_f (median) | temp_bias | |---|---|---|---|---|---| | default | 3 | −0.94 | 2.42 | 5.52 (5) | −0.57 | | pf1.5 | 3 | −2.13 | 1.44 | 5.71 (5) | −1.79 | | pf3.5 | 3 | +0.01 | 3.39 | 5.42 (5) | +0.39 | | mf3 | 3 | +0.91 | 2.40 | 3.25 (3) | +1.30 | | mf7 | 3 | −1.97 | 2.43 | 7.79 (7) | −1.61 | | mf9 | 3 | −2.64 | 2.44 | 9.87 (9) | −2.28 | | centerlines | 11 | −0.97 | 2.43 | 5.50 (5) | −0.59 | | lambda6 | 26 | −0.94 | 2.42 | 5.52 (5) | −0.57 | - With their own temp bias files, the mf presets now hold: the median melt_f equals the preset. - Inversion volume is ~127.9 km³ everywhere (calibrated to consensus). - Extra failures: centerlines 9 `initialize_flowlines` + 2 `glacier_masks`; lambda6 24 small glaciers (1.7 km² in total) in `init_present_time_glacier`, "trapezoid section below its physical minimum". - **Still to do:** - Delete the `exps/*/RGI62/b_160/L3a` scratch (~4 GB each). - Delete the `test_run/` and `_logs_test/` leftovers in `exps/lambda6` and `exps/centerlines`. - **History of today:** - A first round (shared temp bias files, L5) was deleted. The mf builds reused the mf5 temp bias file, so the preset melt_f survived in only ~8-30% of the glaciers: area-weighted melt_f spanned 5.0-6.5, not 3-9. Hence one temp bias run per experiment. - Two later submissions were cancelled by `afterok` when their smoke test failed on a GitHub download (salem/OGGM sample data, once an empty file, once an HTTP 503). - **Sample data:** `_sample_data/` holds the salem (`8ad39d15...`, the container's) and OGGM (`ae1dfa73...`, master's) sample data. `run_prepro.sh` symlinks it into each job's fake `$HOME`, so no job downloads from GitHub. Update it if the OGGM commit or the container changes, otherwise the jobs quietly go back to downloading. - **OGGM PR:** the user is upstreaming the two hunks of `oggm_informed_fixed_prcp_fac.patch`: preset prcp_fac with informed_threestep, and the `get_lock` crash for RGI62/70G. Once merged, set `OGGM_COMMIT` to the merge commit and `OGGM_PATCH=` in `common.sh`. - **Not done (forcing paths):** the forcing file paths are not recorded in the gdirs. They were checked by hand: a test gdir's climate file matches `monthly_hists_v20261001/ensemble/tp_monthly_ens183`, `ensemble/p_monthly_ens183` and `elev.nc` exactly. Writing the paths as attributes in `lmr_forcing.py` was offered. ## round4: historical runs from the v5 gdirs ### Crew meeting 2026-10-19: paper plan and next steps Full notes and tasks are in the user's Firn vault (`Projects/Alps attribution`). - Paper theme: timing of emergence vs total industrial-era attribution. Working title "Little Ice Age glacier state sets when anthropogenic glacier loss emerges in the Alps". - Next steps: - a script for Mira to reproduce the dip in the NAT 30-yr spread: done, `climate_analysis/nat_spread_for_mira/nat_spread_hef.py` (+ png). It reads her stacked `TpNAT_monthly_ens###.nc` directly and uses the HEF cell 10.5°E 46.5°N, JJAS, 31-yr centred mean. The std across members goes ~0.10 K -> ~0.045 K around 1865 -> ~0.22 K around 1935. It runs in ~20 s. **Sent to Mira 2026-10-19** (+ an untested `.m` version); - a summary email to the crew: **sent 2026-10-19**. It asks Mira for HIST strands of all members from 1850 on, or else Tier A / Tier B below; - **Uncertainty sources of the NAT spread (checked 2026-10-19, glacier-weighted JJAS, newnat):** - Pre-1850, the 31-yr spread (0.109 K) ≈ pure white noise: IAV 0.59/√31 = 0.105. So LMR's own low-frequency member spread contributes only ~0.03 K, and the members' 1500–1849 means differ by only 0.034 K (std). There is no decadal persistence (lag-1 ac ~0). - Post-1900, the spread is 0.21 K: member offsets (std of the 1901–2004 means) 0.16 K, plus member trends (std 0.28 K/century), plus a residual of 0.11 K (= white noise). IAV is rescaled to obs and white (lag-1 ~0.06). - So the extra spread is structural (per-member level and trend), not interannual or decadal variability. - All 197 post-1900 series are distinct, with no pair at r > 0.99. The median pairwise r of the annual series is 0.03: no shared volcanic/solar signal is visible, so volcanoes are a minor contributor. - The 1961–1990 means differ by 0.20 K std (11.01–12.31 °C), and are uncorrelated with the member's 1820–1850 mean (r −0.1). The GCM level is independent of the LMR past. - Guess (to ask Mira): the GCM anomaly reference / anchoring, model drift and different responses give per-member offsets that are not natural variability. They inflate the post-1900 SNR noise, so the SNR is conservative. They are likely behind the "crazy warm" members (ens203 +0.8 K in 1961–1990). - a new HIST member list for Mira. **Decided (user, 2026-10-19): ask for both rankings, 9 strands:** 1850 volume quantiles ens066 / 166 / 044 / 197 / 139, and 1820–1850 mean quantiles ens066 / 008 / 010 / 055 / 041. Figure: `round4/attribution/member_choice.py` -> `figs/member_choice_lam8_jjas_dyn_newnat.png`; - one more OGGM calibration round (a few members should come close to the observations); - why some NAT members are very warm and pull the SNR down; - emergence vs total loss vs SNR; - a commitment experiment for all HIST and NAT (committed volume loss); - writing. ### RESUME HERE (2026-10-07) **DECISION (user, 2026-10-07): for the next steps (attribution) we use the JJAS configuration, `lam8_jjas_dyn`.** The monthly one (`lam8_monthly_dyn`) is kept for comparison only. **Attribution dataset = `round4/ensemble/output/lam8_jjas_dyn/`** (see the README.md there): - Per-glacier annual runs, years 1–2025, 3892 glaciers (11-01): - `lmr_mira_hist_ens{053,123,124,153,178,183}_reg11.nc` (6 HIST; ens123 is in the data but left out of the plots, user's choice) - `lmr_mira_nat_ens###_reg11.nc` (197 NAT; 70 end before 2020: 2006 / 2010 / 2013 / 2015, as their forcing files) - `lmr_mira_hist_reg11.nc` -> ens183 (the calibration member) - Regional CSVs on the common subset (3660 glaciers, 99.4% of the area): `round4/ensemble/filtered/lam8_jjas_dyn/` (`volume_km3_*`, `area_km2_*`, `glacier_statistics_reg11.csv` = the subset ids, `lia_max_members.csv`). - **HIST ens### == NAT ens### bit-for-bit before 1850** (checked for 183, 053, 178): HIST − NAT of the same member isolates the 1850+ forcing. After 1850, HIST and NAT differ only in temperature (precip is the same member file). - Model: v5 gdirs `lam8_jjas` (λ8, consensus, prcp_fac 2.5, melt_f prior 5, informed 3-step with its own temp bias file) + per-glacier dynamic melt_f (`round4/sensitivity/output/lam8_jjas_dyn/dyn_calib_reg11.csv`), calibrated on HIST ens183 only, the same for all members. - **Attribution (2026-10-07, for a talk):** `round4/attribution/attribution.py` (Mira's SNR method, signed, noise = NAT-ensemble spread of the change, 15-yr window vs 1820–1850, same calendar window for HIST/NAT, stops 2021). Env vars EXP / MEMBER / WINDOW. Figures + CSVs in `round4/attribution/figs/`; per-glacier NAT cube cached there (`cube_.nc`; building it takes ~20 min with 8 procs, ~1 min/member CPU-bound netCDF decoding). ens183 jjas regional volume SNR: 1975 1.0, 2000 1.9, 2021 3.9 (outside all 127 NAT). - **How the forcing is built** (Mira, email 2026-10-07; paper in prep, do not distribute): LMR anomalies vs 1960–1990, 30-yr low-pass, flat (850–900 mean) for years 1–850 with a blend over 850–900, plus white noise with the observed interannual std; then 1850–1900 linear blend into BE (HIST) or into a CMIP5/6 hist-nat GCM run (NAT); all anchored on the observed 1960–1990 climatology. NAT keeps the GCM's own variability. Checked (`round4/attribution/nat_spread.py`, JJAS, glacier-area weighted over the 20 grid cells with glaciers): NAT is not pinned to obs 1961–1990 (grid mean 17.34 °C vs HIST 17.67); NAT end years 2005/2009/2012/2014/2020 = CMIP experiment lengths; after 1900 the NAT spread doubles (31-yr mean spread 0.11 -> 0.21 K; interannual std 0.58 -> 0.72 K vs obs 0.68; regional volume std 12 -> 20 km³), so the post-1900 noise in the SNR comes from the GCMs. 30-yr JJAS anomaly vs 1800–1850 (`figs/nat_anomaly_*`): HIST +1.7 K in 1991–2020 vs NAT median +0.15 (P95 +0.5); NAT itself is +0.1–0.3 K above 1800–1850. - Forcing: `data/unstacked/lmr_seas_fullsplit_oldBE/`: `jjas_hists_v20261005/ensemble/tp_monthly_ens###.nc` (HIST T), `ensemble/tp_nat_monthly_ens###.nc` (NAT T), `ensemble/p_monthly_ens###.nc` (P, both), `elev.nc`. - **Overnight (submitted 2026-10-05 ~23:20): full rebuild on Mira's fixed HIST files, in two configurations**, each one afterok chain: | config | v5 gdirs (tb1 / build1 / build2) | dyn calib (run / merge) | ensemble (run / merge) | |---|---|---|---| | `lam8_jjas` (HIST_DIR=jjas_hists_v20261005) | 576215 / 576216 / 576217 | 576218 / 576219 | 576220 / 576221 | | `lam8_monthly` (HIST_DIR=monthly_hists_v20261005) | 576222 / 576223 / 576224 | 576225 / 576226 | 576227 / 576228 | - Same setup as lam8_dyn: λ8, consensus, prcp_fac 2.5, melt_f 5 prior, own temp bias run, dyn melt_f from 1000 (maxiter 40, step min 0.3, 0.2 × Hugonnet error). Ensemble: 6 HIST (053 123 124 153 178 183) + 197 NAT. - Outputs: - gdirs `lmr_gdirs_v5_202610/exps/lam8_{jjas,monthly}/` - `round4/sensitivity/output/lam8_{jjas,monthly}_dyn/` - `round4/ensemble/output/lam8_{jjas,monthly}_dyn/` - When done: - check `sacct`; - `cd round4/ensemble && python ../preprocess_mira.py lam8_jjas_dyn` (same for monthly); - plots: set EXP / GDIR_EXP in `plot_talk.py` and `plot_lia_max.py` (delete the lia_max cache); - check that HIST ens183 in the ensemble == sensitivity `_dyn`; - expect `temp_bias/temp_bias.csv` of lam8_jjas == the old lam8 one (static climate unchanged). - **Status 2026-10-06 ~11:30:** all 974 jobs COMPLETED. Checks: - The temp bias file of lam8_jjas == the old lam8 one (max diff 0) ✓. - The gdir climate of lam8_jjas == new `jjas_hists_v20261005` ✓. - **BUG (mine): the 6 HIST members of both ensembles used the OLD `jjas_hists_v20261001` temperatures.** `submit_ens.sh` read HIST_DIR from the command line, then sourced common.sh -> v5 common.sh, which sets `HIST_DIR=jjas_hists_v20261001` as its default and overwrote it. Detected because HIST ens183 in the ensemble != the dyn calibration run from 1850 (Aletsch grew in 1850–1880 like the old run). - NAT members and dyn calibrations are fine: NAT temperature comes from `ensemble/`, and melt_f from the right CSV. - Fixed (`HIST_DIR_ARG` kept across the source; a MEMBERS subset now writes `members_rerun_*.txt`, so members.txt stays complete). Rerun of the 6 HIST members only: jjas run 577191 / merge 577231, monthly 577232 / 577233. The merges overwrite the wrong HIST files. - **DONE 2026-10-06 ~16:30:** HIST reruns merged (6/6 each). - ens183 ensemble == dyn calibration run regionally to ~0.001 km³; ~20 glaciers differ by ≤ 0.08% (jjas) / 0.25% (monthly). That is float64 vs float32 noise (the new HIST files are float64), not a forcing difference. - `preprocess_mira.py` done for both; common subset 3660 (jjas) / 3668 (monthly) glaciers, ~99.3% of the area. - Plots: `EXP=lam8_{c}_dyn GDIR_EXP=lam8_{c} python plot_talk.py` / `plot_lia_max.py` (both now read EXP / GDIR_EXP from the env) -> `talk/talk_7_*`, `talk/lia_max_{volume,area}_*` for both configs. - **Results** (5 HIST: 053 124 153 178 183; ens123 left out as before): | | jjas | monthly | obs | |---|---|---|---| | LIA vol, 1800–1875 per-glacier max: NAT median / HIST | 232 / 222–252 | 242 / 227–261 | 281 ± 43 | | LIA area: NAT median / HIST | 3610 / 3434–3890 | 3736 / 3483–4015 | 4211 ± 215 | | vol 2015 (HIST) | 105–114 | 103–111 | 99.6 | | dV 2000–2020 (HIST, km³ ice) | −39.5 to −53.7 | −38.7 to −52.8 | ~−36.6 | | area 2003 (HIST) | 2218–2347 | 2186–2313 | ~2045 (RGI) | | dA 2003–2015 (HIST) | −15.9 to −18.1% | −15.4 to −17.6% | −13.5% | - With the fix, the HIST members converge quickly after 1850: no more ~85 km³ spread in 1870–1900 (1900: 197–205 jjas, 183–196 monthly). - Monthly gives a ~10 km³ larger NAT LIA and HIST closer to the reference; present day is about the same for both. - Only ens183 (the calibration member) gets close to Hugonnet; the others lose more (different precip and glacier states). - NAT 1850 percentiles of the HIST members: jjas 124 P14, 053 P18, 183 P41, 153 P78, 178 P98; monthly 124 P14, 053 P21, 183 P45, 153 P78, 178 P98. - **Before this:** - `monthly_hists_v20261001` was renamed `jjas_hists_v20261001` (it was the JJAS variant); `HIST_DIR` in the v5 configs (default `jjas_hists_v20261001`), `HIST_DIR` required by `submit_ens.sh`. - All round4 outputs of the v20261001 era moved to `round4/_archive_v20261001/` (sensitivity output + lia_max.csv; ensemble output, filtered, talk). - **Goal:** find the best simulation, then do sensitivity around it. - **Priority:** LIA area and volume close to Reinthaler (1800–1900 per-glacier max: 4211 km² / 281 km³), without sacrificing the 2015 area and Hugonnet 2000–2020 too much. - Centerlines are dropped from now on (a pain with higher λ). - **Done:** per-glacier dynamic melt_f calibration on `lam8_ib` -> `round4/ sensitivity/output/lam8_ib_dyn/` (jobs 574550/574551, all OK; the slowest chunk took 36 min). - **Convergence:** - Full success: 2838 glaciers (67% of the area). - Part success: 777 (29%). Of these, 420 hit maxiter 20 and 116 hit melt_f 17. - Fallback ("no spinup"): 176. First guess only: 59 (3% of the area). - Median |mismatch| went from 681 to 43 kg m⁻² yr⁻¹; 77% are within tolerance. - **melt_f:** median 5.00 -> 5.42; the area-weighted mean is unchanged (5.52 -> 5.51). 63% went up, 26% down. - **Result against lam8_ib:** - Mass change 2000–2020: −33.6 Gt (Hugonnet −33.0), from −38.3. - 2003 area: 2347 km² (+15%), from 2481. - LIA area: 4155 km² (−1%). LIA volume: 231 km³ (−18%, unchanged). - Area ratio 1.77, volume ratio 2.40. - Area change 2003–2015: −15.4% (−13.5% observed). - So the LIA was *not* broken by the calibration. - That output was renamed to `lam8_ib_dyn20` (maxiter 20). The `_dyn` names are now maxiter 40. - **Done** (2026-10-04 ~12:30): dynamic calibration with maxiter 40 on 4 bases (runs 574591/3/5/7, merges 574592/4/6/8, all OK, the slowest chunk 70 min). Against their uncalibrated bases (1800–1900, 3438 glaciers): | exp | LIA area (4211) | LIA vol (281) | 2003 area (2042) | dA 03–15 (−13.5%) | dM 00–20 (−33.0) | ratios (1.99 / 2.82) | |---|---|---|---|---|---|---| | lam8_ib_dyn | 4173 | 234 | 2366 | −15.1% | −33.9 | 1.76 / 2.38 | | lam8_dyn | 4016 | 260 | **2299** | −15.1% | −35.8 | 1.75 / 2.24 | | lam8_ib_mf7_dyn | 4555 | 255 | 2460 | −16.1% | −35.8 | 1.85 / 2.42 | | lam6_ib_dyn | 3980 | 231 | 2336 | −15.2% | −33.8 | 1.70 / 2.40 | - With consensus (lam8), the calibration *does* lower the LIA: 4622 / 282 -> 4016 / 260. - **OGGM bug found** in `run_dynamic_melt_f_calibration` (`core/dynamic_spinup.py:2402`). The max step per iteration uses `[melt_f_min, melt_f_min]` instead of `[melt_f_min, melt_f_max]`. - So the step is max(|1.5 − melt_f_initial| / maxiter, 0.1), and melt_f can only move by ~±3.5 in total from 5, whatever maxiter is. - That is why maxiter 40 did not help (full successes 2838 -> 2840): every part success moved by at most iterations × step cap. ~800 glaciers (~30% of the area) are part successes, ~600 of them still losing too much mass. - **Fixed** as a third hunk in `oggm_informed_fixed_prcp_fac.patch` (backup of the 2-hunk version: `.bak_2hunks`). The user sent a bug report to another Claude for the upstream fix. - Local test, chunk 39: mixed. One glacier newly converged, a couple are slightly worse (larger steps on a non-monotonic response); full successes went 16 -> 15. - The buggy outputs were renamed to `*_dynbug`. **Rerunning** with the fix (maxiter 40, so the step is 0.3) since 2026-10-04 ~13:00: | output | run | merge | |---|---|---| | `lam8_ib_dyn` | 574755 | 574756 | | `lam8_dyn` | 574757 | 574758 | | `lam8_ib_mf7_dyn` | 574759 | 574760 | | `lam6_ib_dyn` | 574761 | 574762 | - When done: `python lia_max.py` (`dyn_calib` table: base / dynbug / dyn), plus the convergence summary of each `dyn_calib_reg11.csv` against its dynbug version. - **Bug report review** (by the OGGM Claude, 2026-10-04): - The bug is real (since #1527), but the fix changes all dynamically calibrated glaciers with melt_f_initial < ~9.25, so it breaks prepro reproducibility. - It also re-enables the error-recovery search (L2667-2673), which today rounds to the same value whenever the step is < 0.5. - Its options: A = fix only; B = fix + an explicit step argument (recommended); C = keep the cap as is and document it. - My earlier evidence was tautological. The proper test (|Δmelt_f| ≈ (n−1)·step and still outside tolerance) gives only 138–176 truly reach-limited glaciers (0.7–1.9% of the area). ~80% of the ~800 part successes fail for other reasons (errors at the limits, repeated guesses, non-monotonic response). - |Δmelt_f| > 3 affects 5–7% of the area. - **My proposal to the user:** - Upstream: B, but with a separate argument for the max deviation from the initial melt_f, its default reproducing today's cap. - Here: drop the 3rd patch hunk and pass `melt_f_max_step_length_minimum` (already a kwarg) from `run_hist.py`; then test the cap (0.1 / 0.3 / 0.6) as a prior. - Optionally, classify the non-reach part successes. - **Done** (user: "whatever you think is right"): - The OGGM patch is back to its 2 hunks. The 3-hunk version is kept as `lmr_gdirs_v5_202610/oggm_patch_with_dynstep_fix.patch.ref`. The running `_dyn` jobs (574755–574762) had all started, i.e. had already applied the 3-hunk patch, before the swap. - From now on: `DYN_STEP_MIN=0.3` in `round4/common.sh` -> `melt_f_max_step_length_minimum`. The dynbug runs ≈ 0.1, the dyn runs ≈ 0.3. - Chunk 39 test, arg 0.3 against the patched formula: same convergence counts, but 6/23 glaciers end on a different melt_f (up to 2.35). The result is path-dependent (non-unique within tolerance). - **Part successes classified** (lam8_ib_dynbug, 807 glaciers): | category | n | % of area | |---|---|---| | maxiter reached, oscillating / non-monotonic | 190 | 12.9 | | stopped early, other (repeated guess / error) | 145 | 8.8 | | stopped early, glacier exceeds domain | 234 | 7.5 | | at melt_f bound | 129 | 0.8 | | reach-limited (step cap) | 109 | 0.5 | - The bug barely matters for us. The levers are the non-monotonic response and the domain exits in the 1000->2025 runs. Border 160 is the largest elev-bands L1/L2 border available. - **Done** (2026-10-04 ~14:15): the 4 `_dyn` reruns with the fixed step (574755–574762, all OK). As predicted, the fix changes little. - Full-success area rises by 0.3–4.7 points; |Δmelt_f| > 3 grows a little. - Regional metrics change by < 1% from the dynbug runs, except lam6_ib: its LIA volume goes 231 -> 244 (path dependence on a few large glaciers?). | exp | LIA area (4211) | LIA vol (281) | 2003 area (2042) | dA 03–15 (−13.5%) | dM 00–20 (−33.0) | ratios (1.99 / 2.82) | |---|---|---|---|---|---|---| | lam8_ib_dyn | 4167 | 235 | 2363 | −15.2% | −33.8 | 1.76 / 2.36 | | lam8_dyn | 3989 | 256 | 2279 | −15.3% | −35.6 | 1.75 / 2.27 | | lam8_ib_mf7_dyn | 4559 | 256 | 2463 | −16.0% | −35.8 | 1.85 / 2.42 | | lam6_ib_dyn | 4007 | 244 | 2360 | −15.2% | −33.8 | 1.70 / 2.25 | - **Full ensemble dataset on lam8_dyn** (user's choice of base, 2026-10-04). Submitted ~15:15: array 574923 (440 tasks = 40 chunks × 11 batches of 20 members, 8 cpus, 8 h), merge 574924. - Output: `round4/ensemble/output/lam8_dyn/`: - `lmr_mira_hist_ens{053,123,124,153,178,183}_reg11.nc` - `lmr_mira_nat_ens###_reg11.nc` (197 members) - `lmr_mira_hist_reg11.nc` -> the ens183 link - `members.txt` - About 100 MB per member, 20 GB in total. - How it works (`submit_ens.sh` -> `run_slurm.sh ... ens` -> `run_ens.py`; merged by `merge_ens.py`, which also deletes the chunks): - The gdirs are v5 `lam8`. The dynamically calibrated melt_f from `sensitivity/output/lam8_dyn/dyn_calib_reg11.csv` is written back into them (static melt_f where there is none). - Then, per member: write its climate (`lmr_forcing.process_mira_data`, HIST temp from `monthly_hists_v20261001`, NAT from `tp_nat_monthly`) and run from year 1. - The calibration stays the ens183 HIST one, as in round3. - Check: HIST ens183 through this pipeline is bit-identical to `sensitivity/output/lam8_dyn` (local and cluster tests on chunk 39). - **DONE (2026-10-04 ~22:00): the full dataset is complete.** - 203/203 members (6 HIST, 197 NAT) in `round4/ensemble/output/lam8_dyn/` (19 GB, no chunks left), plus round3-style CSVs in `round4/ensemble/filtered/lam8_dyn/` (`preprocess_mira.py lam8_dyn`). - Final common subset (hist ∩ nat ∩ hist_ens): 3653 glaciers, 99.2% of the area. The 1850 volume across HIST members is 192–253 km³. - NAT members end in 2006/2010/2013/2015/2021 (8/3/36/23/127 members): that is when their `tp_nat_monthly` files end, identical in round3. A data property, not a bug. - **The merge slowness was not a missing `.load()`.** merge_ens.py started a new python per member (subprocess); on the compute nodes every one re-imported xarray from NFS, which takes minutes. One member merges in ~12 s in-process. merge_ens.py now imports `merge_chunks.merge()` and runs it in a ProcessPool (`MERGE_JOBS`). The last 129 members were merged on the login node in 4 min. - Figure: `round4/ensemble/plot_talk.py` -> `talk/talk_7_lam8_dyn.png` (round3's talk_7 style, without ens123; the members are ranked by the 1850 volume percentile of their NAT counterpart in this dataset: 124 P12, 053 P20, 183 P42, 153 P79, 178 P98). - LIA figures: `round4/ensemble/plot_lia_max.py` -> `talk/lia_max_{volume,area}_lam8_dyn.png`. Left: the time series (the window shaded). Right: the sum of per-glacier maxima over **1800–1975** for NAT (box), HIST and the reference. **The user's window is 1800–1875** (1975 was a misunderstanding; all three windows are cached). - Per-member values are cached in `filtered/lam8_dyn/lia_max_members.csv` (1800–1900 and 1800–1975; delete it to recompute, ~2.5 min). - Over 1800–1875, volume: NAT median 230 (P5–P95 213–248), HIST 228–255, ref 281. Area: NAT median 3611 (3348–3867), HIST 3563–3990, ref 4211. - The reference is now drawn with its 1-sigma uncertainty from the paper (`REF_AREA_ERR=0.051`, `REF_VOL_ERR=0.153`): area 4211 ± 215 -> 3996–4426, volume 281 ± 43 -> 238–324. - Volume: the HIST members are mostly inside the range; only ens124 and ens153 (228–231) are just below. - Area: all members are below the range; ens183 (3990) is at its edge. **The area deficit is the robust model shortfall.** - **Status 21:45:** all 440 run tasks COMPLETED. The first merge (574924) was too slow: ~2.5 min/member serially (40 chunks of ~2.5 MB on NFS), so it would have hit its 4 h limit. Cancelled after 73 of 203 members. - `merge_ens.py` now merges `MERGE_JOBS` (8) members in parallel, and submit_ens.sh gives the merge 8 cpus and 12 h. - The remaining 129 were resubmitted as merge job **575371** (log `round4/ensemble/_logs/lam8_dyn_merge2_575371.out`). - Restart-safe: merged members' chunk folders are deleted, the others are untouched (no `.tmp` was left). - When done: - `cat round4/ensemble/_logs/lam8_dyn_merge*.out` - then `cd round4/ensemble && python ../preprocess_mira.py lam8_dyn` (round3-style CSVs in `ensemble/filtered/lam8_dyn/`). - Recovery if tasks failed: `TASKS= MELT_F_CSV=... ./submit_ens.sh ensemble lam8 lam8_dyn`, then `$LOCAL_PY merge_ens.py ensemble/output/lam8_dyn 40`. Merged members are skipped, because their chunk folders are already gone. - **Best so far:** - lam8_ib_dyn: LIA area −1%, Hugonnet matched, but LIA volume −17%. - lam8_dyn: LIA area −5%, volume −9%, the smallest 2003 overshoot. - Still too large in 2003 everywhere (+12–21%), area loss 2003–2015 too fast (−15 to −16%), area ratio ≤ 1.85. - **Domain exits during the dyn calibration are a symptom, not a border problem** (checked on lam8_ib_dyn; lam8_dyn is the same): - 432 glaciers (8.8% of the area) exit their domain at some point. - 176 fall back and have no output, but they have no output in the uncalibrated base run either: tiny glaciers that leave early. - 256 keep their output: the best melt_f that ran, or the static one. - Glaciers without output total only 0.6% of the area, calibrated or not. - Those 256 are still ~2× their RGI area in 2003 (median; others 1.09×) and keep a large residual mismatch (median −734 kg m⁻² yr⁻¹). They never retreat enough, whatever melt_f is tried. - Their static calibration looks like everyone else's (temp_bias −0.57, prcp_fac at the ±20% bound for ~70%, melt_f 5). - Also: the *median* glacier reaches ~6× its RGI area at its 1000–1900 max. Small glaciers overgrow a lot; the regional area ratio (1.77) is dominated by the big ones. The overgrowth of small glaciers is the thing to understand. - Do NOT enlarge the border for this: it would only let them grow more. - **Ideas not yet done:** - Calibrate melt_f per 1° grid point instead of per glacier, which is more robust to Hugonnet noise on small glaciers. - λ=10. - Write the forcing paths into the gdirs. - Clean up `exps/*/RGI62/b_160/L3a` (~4 GB each) and the test_run leftovers. - **Rules learned the hard way:** - Each config gets its own temp bias file. - Gdirs stop at L3, border 160. - 11-02 (Pyrenees) is excluded. - Always submit from `round4/` with the study name (`./submit.sh sensitivity ...`). Submitting from inside the study folder creates `sensitivity/sensitivity`. Runs (not gdirs) from the v5 L3 gdirs: start from the model flowlines at the RGI date, in year 1 of the forcing, to 2025. The spinup centuries get discarded later. The climate already in the gdirs (ens183 HIST) is used as-is. Chunked: 40 chunks of 100 glaciers per experiment (array job), then a merge job into `/output//lmr_mira_hist_reg11.nc`, which uses round3's name and layout, so `preprocess_mira.py` (copied unchanged) works. Same OGGM + patch as v5 (`common.sh` sources the v5 one). Centerlines runs need `evolution_model = FluxBased`: `run_slurm.sh` sets it, because it is not stored in the gdirs. See [round4/README.md](round4/README.md). - `sensitivity/`: one run per v5 experiment. Submitted 2026-10-03 ~21:00 (run array jobs 573940–573952 and 574173, merges 573941–573953 and 574174). A smoke test on chunk 39 (`sensitivity/test_run/`) matched a local run exactly. - Expected: a few tiny glaciers with very low melt_f grow out of their domain around year ~105 and have no output (5 of 28 in chunk 39), as in round3. - **Done** (all 8 merged, 2026-10-03 ~22:00). `sensitivity/lia_max.py` writes `lia_max.csv`, on the 3454 glaciers that ran in all experiments (97.7% of the RGI area), using `area_min_h`: - **LIA state:** for each glacier, the max of volume and area over **1800–1900**, summed. Compared with Reinthaler (Alps LIA 4211 km² / 281 km³). - **Present-day fit:** - Initial area (year 1) against RGI: it matches (−0.13%). - 2003 area against RGI. - Area change 2003→2015 against Reinthaler 2015 / RGI 11-01 (−13.5%). - Mass change 2000–2020 against Hugonnet per glacier (−33.0 Gt). - **Findings:** - In 2003 the model is too large in every experiment (+9 to +29%). - It loses mass too fast in all experiments but mf3 (−34 to −53 Gt). - With 1800–1875, most glaciers peaked at the very end of the window (57–89% of the area in 1871–1875). With 1800–1900 only 0.2–17% do: the median peak year is 1871–1892, depending on the experiment. - Larger windows (default): 1600–1900 gives the same as 1800–1900 (+0.1%), because there is no 17th-century maximum. Over 1250–1900 the model's absolute max is **~1384**: 80–88% of the area peaks in 1350–1424, with about +7% area and +10% volume over 1800–1900. The reference maps only the last (~1850) maximum, so stick to 1800–1900. Volume at the year of max area differs from the independent max volume by <1%. - Notes on the reference are in Claude's memory (lia-reference-reinthaler-paul-2025). - Parameters that improve the LIA (high melt_f) worsen the present-day mass change. - **11-02 (Pyrenees) is excluded explicitly** in round4: `run_hist.py` drops it from the chunks, and `lia_max.py` and `preprocess_mira.py` count only 11-01. It never had output anyway, because it is outside the LMR grid and the gdirs fail mass balance calibration. The v5 gdir builds still include it, which is harmless (35 glaciers, 3.2 km², all failed). - **Towards the best run** (submitted 2026-10-03 ~23:00). Five more v5 experiments, all elev bands from L1 (`exps/lam8`, `lam6_ib`, `lam8_ib`, `lam8_mf7`, `lam8_ib_mf7`), each with its own temp bias run, build and round4 run (merges 574275–574295, behind a smoke test 574268–574270). - `ib` = inversion calibrated to IceBoost (`INV_VOLUME=iceboost` in config.sh): 108.6 km³ for 11-01, against 127.9 km³ for consensus. - Rationale: the default runs are too large in 2000–2020 (+10–17% area) and their specific MB is too negative (−0.95 against Hugonnet's −0.81 m w.e./yr), so the glaciers respond too slowly. Higher λ and thinner ice should speed up the response. λ=6 was the best single change. - Centerlines are dropped from now on. - **Results** (1800–1900 window, 3438 common glaciers): | exp | LIA area (4211) | LIA vol (281) | area ratio (1.99) | vol ratio (2.82) | 2003 area (2042) | dA 03–15 (−13.5%) | dM 00–20 (−33.0 Gt) | |---|---|---|---|---|---|---|---| | lambda6 | 4386 | 279 | 1.68 | 2.36 | 2611 | −13.8% | −45.0 | | lam6_ib | 4209 | 232 | 1.71 | 2.51 | 2463 | −14.7% | −38.0 | | lam8 | 4622 | 282 | 1.75 | 2.38 | 2634 | −14.3% | −45.2 | | lam8_ib | 4426 | 234 | 1.78 | 2.53 | 2481 | −15.2% | −38.3 | | lam8_mf7 | 5081 | 311 | 1.83 | 2.45 | 2774 | −15.8% | −50.8 | | lam8_ib_mf7 | 4833 | 257 | 1.85 | 2.57 | 2608 | −16.2% | −42.9 | - IceBoost scales all volumes down by ~15%: the 2015 volume is ~93 km³ (Reinthaler 99.6), against ~118 km³ with consensus. Less overshoot in 2003, less mass loss, better ratios, but the absolute LIA volume is −17%. - λ8 raises the LIA area and the ratios. mf7 overshoots the LIA area and the present-day mass loss. - Candidate base for the per-glacier dynamic melt_f calibration (discussed, not started): lam8_ib. - Four merge jobs timed out (30 min) because `merge_chunks.py` sorted lazily opened files. It now loads the chunks first (10 s). The four were merged locally. - **Per-glacier dynamic melt_f calibration** (`DYN_CALIB=1 ./submit.sh sensitivity lam8_ib`, output `sensitivity/output/lam8_ib_dyn/`, submitted 2026-10-04 ~00:15, jobs 574550/574551). - Method: OGGM's `run_dynamic_melt_f_calibration` with the plain `dynamic_melt_f_run` (no dynamic spinup, i.e. no temp_bias loop). Each iteration runs the glacier from the RGI-date flowlines in 1000 to 2025 and matches its 2000–2020 mass balance to Hugonnet (from the gdir's observations.yml) within 0.2 × its error. Only melt_f changes. - Afterwards, the usual year-1 run uses the calibrated melt_f. `dyn_calib_reg11.csv` has the per-glacier diagnostics. - Gotchas: - The calibration needs `mass` in `store_diagnostic_variables` (now added). - Read the diagnostics from `diagnostics.json`: the main process's `gdir.settings` is a stale copy, because the calibration runs in the workers. - Local test on chunk 39: RGI60-11.03901 went melt_f 5.45 → 6.63 (up!) in 7 iterations. With temp_bias fixed, a higher melt_f gives a smaller glacier and less mass loss. The year-1 and 1000 starts give identical 2000–2020 results for it. ## round4/validation: single-glacier validation (2026-10-12) `round4/validation/validate_hef.py` (EXP / GDIR_EXP / MEMBER from the env): Hintereisferner, `lam8_jjas_dyn` HIST ens183 against WGMS mass balance (`gdir.get_ref_mb_data`) and Leclercq length (`gdir.get_ref_length_data`). Nothing is re-run: it reads the ensemble output, and fetches the v5 gdir (into `_wd/`) only for the observations and the fixed geometry MB. CSVs next to the script, figures in `figs/`. - Each MB is calibrated to Hugonnet 2000–2020 on its own geometry: the transient run with the dyn melt_f (4.7), the fixed 2003 geometry with the gdir's static melt_f (5.0, gives −1100.3 = Hugonnet exactly; with 4.7 it would give −928). Both melt_f are passed explicitly, nothing is written to the gdir settings (user's concern: parameter cross-contamination). - MB 1953–2020 (mm w.e./yr, WGMS −653): transient (dV / area, as the OGGM tutorial) −588, bias +65; fixed geometry −456, bias +197. r 0.76 for both. 2000–2020: bias −16 / +7 (the calibration period). 1953–1979: +115 / +342. Decadal means in `figs/hef_mb_decadal_*.png`. - Length change relative to 2003: good from ~1955 (1975: +500 m vs +389), but the model has only +700 m in 1855 against +2858 m observed. HEF never gets long enough at the LIA (other HIST members: +800 to +1300 m). - MB trend 1953–2020 (mm w.e./yr per decade): WGMS −181 ± 32, transient −206, fixed geometry −251. - Centerlines check (`run_hef_centerlines.py`, PYTHONPATH = the v5 folder for lmr_forcing; own `_wd_centerlines/`; ~1 min): the v5 `centerlines` gdir (default λ, static melt_f 5) re-forced with `jjas_hists_v20261005` ens183 (identical to the gdir's climate from 1900, asserted), FluxBased, from year 1. Length vs 2003: +1100 m in 1855, +1500 in 1900 (max), so somewhat longer than elev bands but still far from the observed +2858, and too long in 1950–1990 (+1300 in 1950 vs +1107, +800 in 1975 vs +389). **The flowline representation is not the main cause of the LIA length deficit.** Both runs are ~700–800 m longer in 2003 than the RGI geometry they start from in year 1 (5000 -> 5800 m bands, 6900 -> 7600 m centerlines). Not like-for-like: λ and static/dyn calibration differ too.