# As attr_fraction_conditional.py, but with the committed loss: the HIST # strands continue with their committed-loss runs (round6 commit dataset: # from the 2025 geometry, random climate 2004-2024, 150 yr, 5 seeds), and # the natural runs are held constant at their final state (they are flat # since ~1900). # # committed HIST volume V_hc = mean over 2081-2100 and seeds (shown up to # 2100; the 2125-2175 equilibrium is ~2 km3 lower) # natural volume V_n,ij = mean of NAT run ij's last 5 years # committed anthropogenic loss A_c = V_n - V_hc, total L_c = V(1850) - V_hc # committed share f_c = A_c / L_c # For t >= 2025, panel (b) follows f(t) along the commitment runs (seed mean) # with the natural volumes held constant. # # Panels: (a) regional volume 1820-2175; (b) share of the loss since 1850, # 1950-2175; (c) committed share vs the 1850 volume (dotted: the 2015 share # of attr_fraction_conditional.py), with Reinthaler (V(1850) = model HIST # 2015 + 180 ± 39 km3, as there). # Output: figs/attr_fraction_committed_.png / .csv, suffix = what # COMMIT adds to commit_pf3.5_mf5_lam12 (e.g. _c2014-2024: random climate # 2014-2024 instead of 2004-2024, read from the name). # States without their own commitment run (the c2014-2024 dataset has ens010 # only) use ens010's: the 7 strands are identical at the end (7-strand run). # Env: NAME, COMMIT (the commitment dataset), ENS. import glob import os import re import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt NAME = os.environ.get('NAME', 'perm_pf3.5_mf5_lam12') COMMIT = os.environ.get('COMMIT', 'commit_pf3.5_mf5_lam12') ENS = os.environ.get('ENS', 'ens_pf3.5_mf5_lam12') STATES = ['123', '023', '117', '054', '010', '041', '178'] COLORS = dict(zip(STATES, plt.cm.Reds(np.linspace(0.35, 0.95, len(STATES))))) BLUES = dict(zip(STATES, plt.cm.Blues(np.linspace(0.35, 0.95, len(STATES))))) m = re.search(r'_c(\d{4})-(\d{4})', COMMIT) CLIM = f'{m.group(1)}–{m.group(2)}' if m else '2004–2024' SUFFIX = COMMIT[len('commit_pf3.5_mf5_lam12'):] if COMMIT.startswith('commit_pf3.5_mf5_lam12') else '_' + COMMIT FALLBACK = '010' # the strand whose commitment run is used if a state has none REINTHALER_DV, REINTHALER_ERR = 180., 39. # Ranges: 16-84% (~ ±1 sigma) everywhere, as attr_fraction_conditional.py; at the # Reinthaler volume NAT and Reinthaler are combined by sampling. QLO, QHI = 0.16, 0.84 QTAG = '16–84%' N_MC = 200_000 EQ = (2081, 2100) # the committed volume: mean over these years (2125-2175: 48.0 km3, ~2 lower) T_STOP = 2100 # end of the commitment shown NAT_LAST = 5 # natural volume: mean of each run's last NAT_LAST years FRAC_START, T_HIST = 1950, 2025 v = pd.read_csv(f'../filtered/{NAME}/volume_km3_perm_reg11.csv', index_col=0) pre = pd.read_csv(f'../filtered/{NAME}/pre1850_reg11.csv', index_col=0) v5 = pd.read_csv(f'../../round5/filtered_ens/{ENS}/volume_km3_nat_ens_reg11.csv', index_col=0) for df in (v, pre, v5): df.index = df.index.astype(int) IDS = pd.read_csv(f'../../round5/filtered_ens/{ENS}/glacier_statistics_reg11.csv', index_col=0).index.values # Commitment runs: regional volume on the same subset cache = f'../filtered/{COMMIT}_volume_km3_reg11.csv' if not os.path.exists(cache): cols = {} for f in sorted(glob.glob(f'../output/{COMMIT}/commit*_state*_seed*_reg11.nc')): s, seed = re.search(r'state(\d{3})_seed(\d+)', f).groups() with xr.open_dataset(f) as ds: cols[f'state{s}_seed{seed}'] = ds.volume.sel(rgi_id=IDS).sum('rgi_id').to_series() * 1e-9 pd.DataFrame(cols).rename_axis('year').to_csv(cache) com = pd.read_csv(cache, index_col=0) com.index = com.index.astype(int) res, rows = {}, [] for s in STATES: nat = v[[c for c in v.columns if c.startswith(f'state{s}_nat_')]] hist = v[f'state{s}_hist'] v0 = hist.loc[1850] c = com[[x for x in com.columns if x.startswith(f'state{s}_')]] if c.shape[1]: assert np.allclose(c.loc[T_HIST], hist.loc[T_HIST], rtol=1e-6) else: c = com[[x for x in com.columns if x.startswith(f'state{FALLBACK}_')]] v_hc_seeds = c.loc[EQ[0]:EQ[1]].mean() v_hc = v_hc_seeds.mean() v_n = nat.apply(lambda x: x.dropna().iloc[-NAT_LAST:].mean()) # per NAT run A_c, L_c = v_n - v_hc, v0 - v_hc f_c = A_c / L_c # the share along time: history, then the commitment run (seed mean), NAT constant h_ext = pd.concat([hist.loc[:T_HIST - 1], c.mean(axis=1)]) n_ext = nat.copy().reindex(h_ext.index) n_ext.loc[2021:] = np.nan n_ext = n_ext.fillna(v_n) # held constant n_ext.loc[:2020] = nat.loc[:2020].where(nat.loc[:2020].notnull(), v_n, axis=1) f_t = n_ext.sub(h_ext, axis=0).div(v0 - h_ext, axis=0) res[s] = dict(hist=hist, com=c, v0=v0, v_n=v_n, v_hc=v_hc, q_t=f_t.quantile([QLO, 0.5, QHI], axis=1).T) rows.append(dict(state=f'ens{s}', V1850=v0, V_hist_2015=hist.loc[2015], V_hist_2025=hist.loc[T_HIST], V_hist_committed=v_hc, V_hist_committed_seed_min=v_hc_seeds.min(), V_hist_committed_seed_max=v_hc_seeds.max(), V_nat_final_median=v_n.median(), committed_loss_since_1850=L_c, committed_anthro_loss_median=A_c.median(), share_committed_plo=f_c.quantile(QLO), share_committed_median=f_c.median(), share_committed_phi=f_c.quantile(QHI))) tab = pd.DataFrame(rows).set_index('state') A_med = tab.committed_anthro_loss_median.mean() Vhc = tab.V_hist_committed.mean() A_q = {p: np.mean([(res[s]['v_n'] - res[s]['v_hc']).quantile(p) for s in STATES]) for p in (QLO, QHI)} Vh15 = tab.V_hist_2015.mean() A15 = np.mean([(v[[c for c in v.columns if c.startswith(f'state{s}_nat_')]].loc[2015] - v[f'state{s}_hist'].loc[2015]).median() for s in STATES]) v_rt = Vh15 + REINTHALER_DV f_rt = A_med / (v_rt - Vhc) f_rt_lo, f_rt_hi = A_med / (v_rt + REINTHALER_ERR - Vhc), A_med / (v_rt - REINTHALER_ERR - Vhc) print(f'committed HIST volume {Vhc:.1f} km3 (2025: {tab.V_hist_2025.mean():.1f}); natural final ' f'{tab.V_nat_final_median.mean():.1f}; committed anthropogenic loss A_c = {A_med:.1f} ' f'(NAT {QTAG} {A_q[QLO]:.1f}-{A_q[QHI]:.1f})') A_all = np.concatenate([(res[s]['v_n'] - res[s]['v_hc']).dropna().values for s in STATES]) rng = np.random.default_rng(0) f_mc = rng.choice(A_all, N_MC) / (rng.normal(v_rt, REINTHALER_ERR, N_MC) - Vhc) f_comb = {p: np.quantile(f_mc, p) for p in (0.05, QLO, 0.5, QHI, 0.95)} print(f'Reinthaler, NAT and Reinthaler combined: median {100 * f_comb[0.5]:.0f}%, {QTAG} ' f'{100 * f_comb[QLO]:.0f}-{100 * f_comb[QHI]:.0f}%, 5-95% {100 * f_comb[0.05]:.0f}-' f'{100 * f_comb[0.95]:.0f}%') print(f'Reinthaler V(1850) = {v_rt:.0f} ± {REINTHALER_ERR:.0f} -> committed share {100 * f_rt:.0f}% ' f'({100 * f_rt_lo:.0f}-{100 * f_rt_hi:.0f}%; NAT {QTAG} at the central value ' f'{100 * A_q[QLO] / (v_rt - Vhc):.0f}-{100 * A_q[QHI] / (v_rt - Vhc):.0f}%); 2015 share: ' f'{100 * A15 / REINTHALER_DV:.0f}%') fig = plt.figure(figsize=(15, 5.2)) gs = fig.add_gridspec(1, 3, width_ratios=[1.25, 1, 1], wspace=0.28) ax_a, ax_b, ax_c = (fig.add_subplot(gs[0, i]) for i in range(3)) t_end = T_STOP com = com.loc[:T_STOP] # (a) volumes for s in STATES: d = res[s] h = pd.concat([pre[f'volume_km3_state{s}'].loc[1820:1849], d['hist'].loc[1850:T_HIST]]) nmed = v[[c for c in v.columns if c.startswith(f'state{s}_nat_')]].median(axis=1).loc[1850:2020] ax_a.plot(h.loc[:1850].index, h.loc[:1850], color='0.5', lw=1) ax_a.plot(nmed.index, nmed, color=BLUES[s], lw=1.6) ax_a.plot([2020, t_end], [d['v_n'].median()] * 2, color=BLUES[s], lw=1.6, ls='--') ax_a.plot(h.loc[1850:].index, h.loc[1850:], color=COLORS[s], lw=1.8) ax_a.plot(d['com'].index, d['com'].mean(axis=1), color=COLORS[s], lw=1.8, ls='--') c10 = res['010']['com'] ax_a.fill_between(c10.index, c10.min(axis=1), c10.max(axis=1), color=COLORS['010'], alpha=0.2, lw=0, label='committed, range of the 5 seeds') ax_a.axvspan(T_HIST, t_end, color='0.5', alpha=0.07, lw=0) ax_a.text(T_HIST + 3, 8, f'committed:\nrandom climate\n{CLIM}', fontsize=7, color='0.35', va='bottom') ax_a.errorbar([1850], [v_rt], yerr=[REINTHALER_ERR], fmt='D', color='k', ms=6, capsize=4, zorder=5, label=f'Reinthaler & Paul (2025), anchored\non model 2015: ' f'{v_rt:.0f} ± {REINTHALER_ERR:.0f} km$^3$') ax_a.plot([], [], color=COLORS['054'], lw=1.8, label='HIST, 7 states (dashed: committed)') ax_a.plot([], [], color=BLUES['117'], lw=1.6, label='NAT from the same 1850 state, median\n(dashed: held constant)') ax_a.plot([], [], color='0.5', lw=1, label='before 1850 (shared)') ax_a.set_xlim(1820, t_end) ax_a.set_ylim(0, 340) ax_a.set_ylabel('Regional volume (km$^3$)') ax_a.set_title('(a) Regional volume: HIST (red) and natural (blue)', loc='left', fontsize=10.5) ax_a.legend(loc='upper right', fontsize=7, frameon=False) ax_a.grid(alpha=0.3) # (b) share since 1850 over time, into the commitment for s in STATES: q = res[s]['q_t'].loc[FRAC_START:] if s == '010': ax_b.fill_between(q.index, 100 * q[QLO], 100 * q[QHI], color=COLORS[s], alpha=0.12, lw=0, label=f'ens010, NAT {QTAG}') ax_b.plot(q.loc[:T_HIST].index, 100 * q.loc[:T_HIST, 0.5], color=COLORS[s], lw=1.8, label=f'ens{s}') ax_b.plot(q.loc[T_HIST:].index, 100 * q.loc[T_HIST:, 0.5], color=COLORS[s], lw=1.8, ls='--') ax_b.axvspan(T_HIST, t_end, color='0.5', alpha=0.07, lw=0) ax_b.axhline(100, color='0.3', lw=0.9) ax_b.plot([], [], color='0.4', lw=1.8, ls='--', label='committed') ax_b.set_xlim(FRAC_START, t_end) ax_b.set_ylim(0, 200) ax_b.set_ylabel('Anthropogenic share of the\nvolume loss since 1850 (%)') ax_b.set_title('(b) Share of the loss since 1850', loc='left', fontsize=10.5) ax_b.legend(loc='upper right', fontsize=7, frameon=False, ncol=2) ax_b.grid(alpha=0.3) # (c) committed share vs the 1850 volume x = np.linspace(Vh15 + 60, 340, 300) ax_c.fill_between(x, 100 * A_q[QLO] / (x - Vhc), 100 * A_q[QHI] / (x - Vhc), color='0.6', alpha=0.2, lw=0, label=f'committed, NAT {QTAG}') ax_c.plot(x, 100 * A_med / (x - Vhc), color='0.3', lw=1.5, label=f'committed: {A_med:.0f} / (V$_{{1850}}$ − {Vhc:.0f})') ax_c.plot(x, 100 * A15 / (x - Vh15), color='0.3', lw=1.2, ls=':', label=f'2015: {A15:.0f} / (V$_{{1850}}$ − {Vh15:.0f})') for s in STATES: d = res[s] f_c = (d['v_n'] - d['v_hc']) / (d['v0'] - d['v_hc']) qlo, q50, qhi = f_c.quantile([QLO, 0.5, QHI]) ax_c.errorbar([d['v0']], [100 * q50], yerr=[[100 * (q50 - qlo)], [100 * (qhi - q50)]], fmt='o', color=COLORS[s], ms=6, capsize=3) ax_c.axvspan(v_rt - REINTHALER_ERR, v_rt + REINTHALER_ERR, color='k', alpha=0.08, lw=0) ax_c.axvline(v_rt, color='k', lw=1, ls='--', label=f'Reinthaler: {v_rt:.0f} ± {REINTHALER_ERR:.0f} km$^3$') ax_c.errorbar([v_rt], [100 * f_comb[0.5]], yerr=[[100 * (f_comb[0.5] - f_comb[0.05])], [100 * (f_comb[0.95] - f_comb[0.5])]], fmt='none', ecolor='k', elinewidth=1, capsize=3, zorder=6) ax_c.errorbar([v_rt], [100 * f_comb[0.5]], yerr=[[100 * (f_comb[0.5] - f_comb[QLO])], [100 * (f_comb[QHI] - f_comb[0.5])]], fmt='D', color='k', ms=6, elinewidth=3, capsize=0, zorder=7, label=f'Reinthaler + NAT: {100 * f_comb[0.5]:.0f}% ({100 * f_comb[QLO]:.0f}–' f'{100 * f_comb[QHI]:.0f}%; 5–95%: {100 * f_comb[0.05]:.0f}–{100 * f_comb[0.95]:.0f}%)') ax_c.axhline(100, color='0.3', lw=0.9) ax_r = ax_c.inset_axes([0, 0, 1, 0.12], sharex=ax_c) ax_r.hist(v5.loc[1850], bins=np.arange(170, 345, 5), color=BLUES['117'], alpha=0.6) ax_r.axis('off') ax_c.text(0.02, 0.135, '1850 volume of the 197 model members', transform=ax_c.transAxes, fontsize=7.5, color=BLUES['054']) ax_c.set_xlim(Vh15 + 60, 340) ax_c.set_ylim(0, 200) ax_c.set_xlabel('Regional volume in 1850 (km$^3$)') ax_c.set_ylabel('Anthropogenic share of the\ncommitted loss since 1850 (%)') ax_c.set_title('(c) Committed share vs the 1850 state', loc='left', fontsize=10.5) ax_c.legend(loc='upper right', fontsize=7.5, frameon=True, framealpha=0.9, edgecolor='none') ax_c.grid(alpha=0.3) fig.suptitle('How much of the Alpine ice loss since 1850 is anthropogenic, including the ' f'committed loss ({CLIM} climate)? NAT conditioned on the 1850 state ({NAME})', y=1.0) fig.savefig(f'figs/attr_fraction_committed_{NAME}{SUFFIX}.png', dpi=150, bbox_inches='tight') tab.round(3).to_csv(f'figs/attr_fraction_committed_{NAME}{SUFFIX}.csv') pd.Series({f'q{int(100 * k):02d}': v for k, v in f_comb.items()}, name='share').round(4).to_csv( f'figs/attr_fraction_committed_{NAME}{SUFFIX}_reinthaler_combined.csv') with pd.option_context('display.width', 250): print(tab.round(2).to_string())