# %% [markdown] # # Iceland Glacier Projections — Comparison with Reference Models # # This notebook compares our Iceland-specific OGGM projections against three # independent reference products, all used in Zekollari et al. (2024, # https://doi.org/10.5194/egusphere-2024-1013): # # | Model | Reference | Notes | # |---|---|---| # | OGGM v1.6.1 standard | Zekollari et al. 2024 | Global default calibration, RGI 6.0 | # | GloGEM | Zekollari et al. 2024 | Independent model, CMIP6 only | # | PyGEM-OGGM | Rounce et al. 2023 | Independent model | # # **Our models** use a custom Iceland-specific calibration and are run for two # historical climate forcings, both using the same mass-balance model (OGGM's default # surface-type-tracking temperature-index model — see `00_experiment_selection_and_ # validation` for how these were chosen): # - **ERA5** (1991–2025), trapezoidal shape, lambda=8 # - **CARRA** (1991–2025), trapezoidal shape, lambda=14 # # **Normalisation:** each model's volume is expressed as a percentage of its own 2020 # value, so that differences in absolute ice volume due to inventory date or calibration # do not affect the comparison — only the *trajectory* (rate of loss) matters. # # **GCM ensemble:** to ensure an apples-to-apples comparison, only the GCMs common # to all models (ours ERA5 ∩ ours CARRA ∩ OGGM std ∩ GloGEM ∩ PyGEM-OGGM) are used. # The pre-2026 segment of our models is the single reanalysis-forced historical run. # # **On data availability:** this notebook ships with its own copy of the data it needs # under `data/` (historical runs, CMIP6 projections, and the reference-model CSVs), so it # can be re-run standalone — it does not depend on the OGGM cluster. # %% import os import yaml import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt # %% [markdown] # ## Configuration # %% data_dir = "data" fig_dir = "figures" os.makedirs(fig_dir, exist_ok=True) our_configs = { "ERA5": { "data_dir": f"{data_dir}/projections/projections_output_era5", "hist_path": (f"{data_dir}/historical/working_dir_adapted_era5_8/" "run_output_spinup_historical.nc"), "label": "OGGM Iceland (ERA5)", "color": "royalblue", }, "CARRA": { "data_dir": f"{data_dir}/projections/projections_output_carra", "hist_path": (f"{data_dir}/historical/working_dir_adapted_carra_14/" "run_output_spinup_historical.nc"), "label": "OGGM Iceland (CARRA)", "color": "firebrick", }, } ref_model_dir = f"{data_dir}/supporting/reference_models" base_oggm_std = f"{ref_model_dir}/oggm_standard" base_pygem = f"{ref_model_dir}/pygem_oggm" base_glogem = f"{ref_model_dir}/glogem" # %% [markdown] # ## Load data # %% for cfg in our_configs.values(): ds_cmip6 = xr.open_dataset(f"{cfg['data_dir']}/merged_2100_projections_cmip6.nc") cfg["vol_cmip6"] = ds_cmip6["volume"].sum(dim="rgi_id") / 1e9 # km³ ds_hist = xr.open_dataset(cfg["hist_path"]) cfg["vol_hist"] = ds_hist.volume.sum(dim="rgi_id") / 1e9 # km³ cfg["mass_hist"] = ds_hist["mass_kg"].sum("rgi_id") / 1e12 # Gt cfg["vol_2020"] = float(cfg["vol_hist"].sel(time=2020)) cfg["ours_norm"] = cfg["vol_cmip6"] / cfg["vol_2020"] * 100 # % of 2020 # %% [markdown] # ## Load and normalise reference models # # Each GCM column is normalised independently to its own 2020 value. # %% def load_ref_norm(base_path, scenario, unit_factor=1.0, start=2020, end=2100): """Load reference CSV (years × GCMs), convert units, normalise each GCM column to its own 2020 value. Returns DataFrame of % relative to 2020.""" df = pd.read_csv(f"{base_path}/{scenario}.csv", index_col=0) * unit_factor df.index = df.index.astype(int) df.columns = df.columns.str.upper() df = df.loc[start:end] return df.div(df.loc[2020]) * 100 ref_sources = { "OGGM standard": (base_oggm_std, 1 / 1e9), # m³ → km³ "GloGEM": (base_glogem, 1.0), # already km³ "PyGEM-OGGM": (base_pygem, 1 / 1e9), } comp_scenarios = ["ssp126", "ssp245", "ssp370", "ssp585"] ref_data = { model_name: {scen: load_ref_norm(base, scen, uf) for scen in comp_scenarios} for model_name, (base, uf) in ref_sources.items() } # %% [markdown] # ## Restrict to common GCMs # # GloGEM and PyGEM-OGGM ran fewer GCMs than our runs and OGGM standard. # All comparisons are restricted to GCMs present in every model (ERA5, CARRA, and all refs). # %% common_gcms = {} for scen in comp_scenarios: gcm_sets = [] for cfg in our_configs.values(): gcms = set(s.split(f"_{scen}")[0].upper() for s in cfg["vol_cmip6"].gcm_scenario.values if s.endswith(f"_{scen}")) gcm_sets.append(gcms) gcm_sets += [set(ref_data[m][scen].columns) for m in ref_sources] common_gcms[scen] = sorted(set.intersection(*gcm_sets)) for model_name in ref_sources: for scen in comp_scenarios: ref_data[model_name][scen] = ref_data[model_name][scen][common_gcms[scen]] def ours_common(cfg, scen): """Return cfg['ours_norm'] subset for scen, restricted to common GCMs.""" idx = [i for i, s in enumerate(cfg["ours_norm"].gcm_scenario.values) if s.endswith(f"_{scen}") and s.split(f"_{scen}")[0].upper() in common_gcms[scen]] return cfg["ours_norm"].isel(gcm_scenario=idx) # %% [markdown] # ## Comparison plot # %% ref_styles = { "OGGM standard": dict(color="seagreen", ls="--", lw=1.8), "GloGEM": dict(color="darkorchid", ls="-.", lw=1.8), "PyGEM-OGGM": dict(color="darkorange", ls=":", lw=1.8), } fig, axes = plt.subplots(2, 2, figsize=(11, 8), sharey=True, sharex=True) axes_flat = axes.flatten() for ax, scen in zip(axes_flat, comp_scenarios): # historical period for each config for cfg in our_configs.values(): hist_years = cfg["vol_hist"].sel(time=slice(2020, 2026)).time.values hist_norm = (cfg["vol_hist"].sel(time=slice(2020, 2026)) / cfg["vol_2020"] * 100).values ax.plot(hist_years, hist_norm, color=cfg["color"], lw=2.5, zorder=5) # our OGGM Iceland projections for each config for cfg in our_configs.values(): ours_scen = ours_common(cfg, scen) med = ours_scen.quantile(0.50, dim="gcm_scenario") p10 = ours_scen.quantile(0.10, dim="gcm_scenario") p90 = ours_scen.quantile(0.90, dim="gcm_scenario") n = ours_scen.sizes["gcm_scenario"] ax.plot(cfg["ours_norm"].time, med, label=f"{cfg['label']} (n={n})", color=cfg["color"], ls="-", lw=2.5, zorder=5) ax.fill_between(cfg["ours_norm"].time, p10, p90, color=cfg["color"], alpha=0.12, zorder=4) # reference models for model_name, style in ref_styles.items(): df = ref_data[model_name][scen] n = len(df.columns) med_r = df.median(axis=1) p10_r = df.quantile(0.10, axis=1) p90_r = df.quantile(0.90, axis=1) ax.plot(df.index, med_r, label=f"{model_name} (n={n})", **style, zorder=3) ax.fill_between(df.index, p10_r, p90_r, color=style["color"], alpha=0.10, zorder=2) ax.axhline(100, color="gray", lw=0.7, ls=":") ax.set_title(scen, fontsize=12) ax.grid(True, alpha=0.3) for ax in axes[1]: ax.set_xlabel("Year", fontsize=11) for ax in axes[:, 0]: ax.set_ylabel("Volume (% of 2020)", fontsize=11) axes_flat[0].legend(fontsize=10, loc="lower left") fig.suptitle("Iceland glacier volume — comparison with reference models (CMIP6)", fontsize=14) plt.tight_layout() fig.savefig(f"{fig_dir}/volume_model_comparison_cmip6.png", dpi=150, bbox_inches="tight") plt.show() # %% [markdown] # ## Mass loss comparison: models vs Hugonnet observations # # Mass change in Gt w.e. for 2000-2020 (the Hugonnet calibration period). # Hugonnet yaml values are in m³ ice (converted from geodetic MB using 850 kg/m³). # Model ice volumes are converted to Gt w.e. using 900 kg/m³ (OGGM default). # For reference models: 10-90 % range across 4 scenarios × common GCMs. # Our models are single reanalysis historical runs (no spread). # GloGEM is excluded here: it only starts in 2015, so it has no 2000-2020 value # to show (see the main comparison plot and the 2100 summary table for GloGEM). # %% with open(f"{data_dir}/supporting/area_dV_obs.yaml") as _f: _obs = yaml.safe_load(_f) RHO_MODEL = 900 / 1000 # km³ ice → Gt w.e. (Models uses 900 kg/m³) def ref_dv_period(y0, y1): """Mass change y0→y1 (Gt w.e.) pooled across all scenarios × common GCMs.""" out = {} for model_name, (base, uf) in ref_sources.items(): vals = [] for scen in comp_scenarios: df_raw = pd.read_csv(f"{base}/{scen}.csv", index_col=0) * uf df_raw.index = df_raw.index.astype(int) df_raw.columns = df_raw.columns.str.upper() if y0 not in df_raw.index: continue for gcm in common_gcms[scen]: if gcm in df_raw.columns: vals.append(float(df_raw.loc[y1, gcm] - df_raw.loc[y0, gcm])) out[model_name] = np.array(vals) * RHO_MODEL if vals else None return out ref_dv_2000_2020 = ref_dv_period(2000, 2020) hug_2000_2020 = _obs["hugonnet_adapted_2000_2020_kg"] / 1e12 hug_2000_2020_err = _obs["hugonnet_2000-01-01_2020-01-01_kg_err"] / 1e12 # mass_kg diagnostic accounts for mixed ice/firn densities — more accurate than volume*0.9 for cfg in our_configs.values(): mh = cfg["mass_hist"] cfg["dv_gt_2000_2020"] = float(mh.sel(time=2020) - mh.sel(time=2000)) # %% # GloGEM omitted from these bar charts — it has no 2000-2020 value (see markdown above). bar_ref_models = ["OGGM standard", "PyGEM-OGGM"] model_names = (["Hugonnet (obs)"] + [cfg["label"] for cfg in our_configs.values()] + bar_ref_models) colors_bar = (["dimgray"] + [cfg["color"] for cfg in our_configs.values()] + ["seagreen", "darkorange"]) def stats(arr): if arr is None or len(arr) == 0: return np.nan, np.nan, np.nan return float(np.median(arr)), float(np.percentile(arr, 10)), float(np.percentile(arr, 90)) our_gt_2000 = [cfg["dv_gt_2000_2020"] for cfg in our_configs.values()] meds = [hug_2000_2020, *our_gt_2000, *[stats(ref_dv_2000_2020[m])[0] for m in bar_ref_models]] p10s = [hug_2000_2020 - hug_2000_2020_err, *our_gt_2000, *[stats(ref_dv_2000_2020[m])[1] for m in bar_ref_models]] p90s = [hug_2000_2020 + hug_2000_2020_err, *our_gt_2000, *[stats(ref_dv_2000_2020[m])[2] for m in bar_ref_models]] fig, ax = plt.subplots(figsize=(7, 5)) y = np.arange(len(model_names)) meds = np.array(meds, dtype=float) p10s = np.array(p10s, dtype=float) p90s = np.array(p90s, dtype=float) xerr_lo = np.where(np.isnan(meds), 0, meds - p10s) xerr_hi = np.where(np.isnan(meds), 0, p90s - meds) ax.barh(y, np.where(np.isnan(meds), 0, meds), xerr=[xerr_lo, xerr_hi], color=colors_bar, alpha=0.8, height=0.6, error_kw=dict(ecolor="black", lw=1.5, capsize=4)) for i, v in enumerate(meds): if np.isnan(v): ax.text(0, i, " n/a", va="center", fontsize=9, color="gray") ax.axvline(0, color="black", lw=0.8) ax.set_yticks(y) ax.set_yticklabels(model_names) ax.set_xlabel("Mass change (Gt w.e.)") ax.set_title("Mass loss 2000-2020") ax.grid(True, axis="x", alpha=0.3) fig.suptitle("Iceland glacier mass loss: models vs Hugonnet observations\n" "(ref models: 10-90 % range across 4 scenarios × common GCMs)", fontsize=12) plt.tight_layout() fig.savefig(f"{fig_dir}/mass_loss_comparison.png", dpi=150, bbox_inches="tight") plt.show() # --- printed table (Gt) --- print(f"\n{'Model':<30} {'2000-2020':>10} {'':>8} {'':>8} (Gt w.e.)") print(f"{'':30} {'median':>10} {'p10':>8} {'p90':>8}") print("-" * 60) for i, name in enumerate(model_names): m1, l1, h1 = meds[i], p10s[i], p90s[i] def fmt(v): return f"{v:10.1f}" if not np.isnan(v) else f"{'n/a':>10}" print(f"{name:<30} {fmt(m1)} {fmt(l1)} {fmt(h1)}") # %% [markdown] # ### Same comparison in km³ ice (no density conversion) # # Uses model ice volumes directly and Hugonnet's original `_m3` values (not adapted). # Allows direct comparison without any density assumption. # %% def ref_dv_period_km3(y0, y1): """Volume change y0→y1 (km³ ice) pooled across all scenarios × common GCMs.""" out = {} for model_name, (base, uf) in ref_sources.items(): vals = [] for scen in comp_scenarios: df_raw = pd.read_csv(f"{base}/{scen}.csv", index_col=0) * uf df_raw.index = df_raw.index.astype(int) df_raw.columns = df_raw.columns.str.upper() if y0 not in df_raw.index: continue for gcm in common_gcms[scen]: if gcm in df_raw.columns: vals.append(float(df_raw.loc[y1, gcm] - df_raw.loc[y0, gcm])) out[model_name] = np.array(vals) if vals else None return out ref_km3_2000_2020 = ref_dv_period_km3(2000, 2020) hug_km3_2000_2020 = _obs["hugonnet_2000-01-01_2020-01-01_m3"] / 1e9 hug_km3_2000_2020_err = _obs["hugonnet_2000-01-01_2020-01-01_m3_err"] / 1e9 for cfg in our_configs.values(): vh = cfg["vol_hist"] cfg["dv_km3_2000_2020"] = float(vh.sel(time=2020) - vh.sel(time=2000)) # %% our_km3_2000 = [cfg["dv_km3_2000_2020"] for cfg in our_configs.values()] meds_km3 = [hug_km3_2000_2020, *our_km3_2000, *[stats(ref_km3_2000_2020[m])[0] for m in bar_ref_models]] p10s_km3 = [hug_km3_2000_2020 - hug_km3_2000_2020_err, *our_km3_2000, *[stats(ref_km3_2000_2020[m])[1] for m in bar_ref_models]] p90s_km3 = [hug_km3_2000_2020 + hug_km3_2000_2020_err, *our_km3_2000, *[stats(ref_km3_2000_2020[m])[2] for m in bar_ref_models]] fig, ax = plt.subplots(figsize=(7, 5)) meds_km3 = np.array(meds_km3, dtype=float) p10s_km3 = np.array(p10s_km3, dtype=float) p90s_km3 = np.array(p90s_km3, dtype=float) xerr_lo = np.where(np.isnan(meds_km3), 0, meds_km3 - p10s_km3) xerr_hi = np.where(np.isnan(meds_km3), 0, p90s_km3 - meds_km3) ax.barh(y, np.where(np.isnan(meds_km3), 0, meds_km3), xerr=[xerr_lo, xerr_hi], color=colors_bar, alpha=0.8, height=0.6, error_kw=dict(ecolor="black", lw=1.5, capsize=4)) for i, v in enumerate(meds_km3): if np.isnan(v): ax.text(0, i, " n/a", va="center", fontsize=9, color="gray") ax.axvline(0, color="black", lw=0.8) ax.set_yticks(y) ax.set_yticklabels(model_names) ax.set_xlabel("Volume change (km³ ice)") ax.set_title("Volume loss 2000-2020") ax.grid(True, axis="x", alpha=0.3) fig.suptitle("Iceland glacier volume loss: models vs Hugonnet observations\n" "(ref models: 10-90 % range across 4 scenarios × common GCMs)", fontsize=12) plt.tight_layout() fig.savefig(f"{fig_dir}/volume_loss_comparison.png", dpi=150, bbox_inches="tight") plt.show() # --- printed table (km³) --- print(f"\n{'Model':<30} {'2000-2020':>10} {'':>8} {'':>8} (km³ ice)") print(f"{'':30} {'median':>10} {'p10':>8} {'p90':>8}") print("-" * 60) for i, name in enumerate(model_names): m1, l1, h1 = meds_km3[i], p10s_km3[i], p90s_km3[i] print(f"{name:<30} {fmt(m1)} {fmt(l1)} {fmt(h1)}") # %% [markdown] # **Key observations:** # - Our Iceland-specific OGGM calibration projects systematically *more mass loss* than # the three reference models across all scenarios, most pronounced under ssp585. # - OGGM standard (same glacier model, global default calibration) is the closest # comparator; the offset reflects the Iceland-specific calibration and the # use of the 2025 inventory with updated bed maps. # - GloGEM consistently retains the most ice, while PyGEM-OGGM sits between OGGM standard # and OGGM iceland. # - The 10-90% spread across the common GCMs is broadly similar across models, # suggesting the inter-model spread is dominated by the GCM forcing rather than # model differences. # %% [markdown] # ## Volume at 2100 — summary table # # Absolute volumes in km³ (median, 10th-90th percentile) across the common GCMs. # %% rows = [] for scen in comp_scenarios: for cfg in our_configs.values(): s = ours_common(cfg, scen).sel(time=2100).values * cfg["vol_2020"] / 100 rows.append(dict(scenario=scen, model=cfg["label"], n=len(s), median=round(float(np.median(s)), 1), p10=round(float(np.percentile(s, 10)), 1), p90=round(float(np.percentile(s, 90)), 1))) for model_name, (base, uf) in ref_sources.items(): df_raw = pd.read_csv(f"{base}/{scen}.csv", index_col=0) * uf df_raw.index = df_raw.index.astype(int) df_raw.columns = df_raw.columns.str.upper() s = df_raw.loc[2100, common_gcms[scen]] rows.append(dict(scenario=scen, model=model_name, n=len(s), median=round(s.median(), 1), p10=round(s.quantile(0.10), 1), p90=round(s.quantile(0.90), 1))) df_vol = pd.DataFrame(rows).set_index(["scenario", "model"]) print(df_vol.to_string()) # %% [markdown] # ## Common GCMs # # GCMs present in ERA5 ∩ CARRA ∩ OGGM standard ∩ GloGEM ∩ PyGEM-OGGM. # %% for scen in comp_scenarios: print(f"{scen}: {common_gcms[scen]}")