# %% [markdown] # # Water divides vs ice divides: Iceland glacier projections # # This notebook compares two OGGM runs that differ **only in how the ice caps are # divided into glacier entities**: # # - **Ice divides (RGI2025)**: the inventory of the first delivery (July 2026), # 572 entities following ice-flow basins. # - **Water divides (RGI2025-hydro v2)**: 571 entities. The ice caps with # Landsvirkjun water divides (Vatnajökull, Hofsjökull, Langjökull, Mýrdalsjökull, # Drangajökull) are split into 59 hydrological catchments, plus 20 "undelineated" # remnants that are not covered by any catchment. All other glaciers are identical # to RGI2025. The total glacier area is the same in both inventories. # # Everything else is identical: CARRA climate, trapezoidal bed (lambda=14), the # regional parameters of the first delivery (precipitation factor, temperature bias, # Glen A), the same OGGM versions and the same CMIP6 runs. Only the melt factor (per # entity), the ice thickness inversion and the dynamic spinup are recomputed for the # new entities. Any difference below is therefore caused by the change of divides. # # We compare CMIP6 projections to 2100. Lines are the median across GCMs, shaded # areas the 10th-90th percentile range. # # This notebook ships with its own copy of the data it needs, under `data/`. # %% import json import matplotlib.pyplot as plt import numpy as np import pandas as pd import xarray as xr # %% [markdown] # ## Load data # %% data_dir = "data" fig_dir = "figures" hist = { "ice": xr.open_dataset(f"{data_dir}/historical/rgi2025/run_output_spinup_historical.nc"), "water": xr.open_dataset(f"{data_dir}/historical/hydro_v2/run_output_spinup_historical.nc"), } proj = { "ice": xr.open_dataset(f"{data_dir}/projections/rgi2025/merged_2100_projections_cmip6.nc"), "water": xr.open_dataset(f"{data_dir}/projections/hydro_v2/merged_2100_projections_cmip6.nc"), } with open(f"{data_dir}/inventory/icecap_ids_rgi2025.json") as f: icecap_ids = {"ice": json.load(f)} with open(f"{data_dir}/inventory/icecap_ids_rgi2025hydro.json") as f: icecap_ids["water"] = json.load(f) labels = {"ice": "Ice divides (RGI2025)", "water": "Water divides (hydro)"} # %% [markdown] # A few low-emission runs that continue to 2300 lose many small glaciers when these # regrow beyond their model domain (a standard OGGM limitation, present in both # inventories). These glaciers are then missing for the whole run, including the # 2026-2100 period. For a clean comparison we exclude these four runs from both # inventories. # %% exclude = ["CESM2-WACCM_ssp126", "CESM2-WACCM_ssp534-over", "MRI-ESM2-0_ssp126", "MRI-ESM2-0_ssp534-over"] runs = [g for g in proj["ice"].gcm_scenario.values if g not in exclude and g in proj["water"].gcm_scenario.values] proj = {k: ds.sel(gcm_scenario=runs) for k, ds in proj.items()} print(f"{len(runs)} CMIP6 runs used in both inventories") # %% [markdown] # ## Why the ice volume differs between the two runs # # Both inventories cover exactly the same glacier area, yet their ice volumes differ # (by about 2% for Iceland, more for individual ice caps, see below). This is # expected: # # - OGGM does not observe ice thickness: it computes it with a mass-conservation # inversion. Along the flowline of each entity, the ice must carry the mass # gained upstream (accumulation) down to where it is lost (ablation). The # thickness needed for that depends on this mass flux, and on the slope and width # of the entity, through the ice creep parameter (Glen A). # - Dividing an ice cap differently changes all of these. A long entity that # collects ice from a large accumulation area carries a larger flux and comes out # thicker; cutting the same ice into smaller pieces gives thinner ice. In reality, # ice flows across water divides, which a catchment-based entity cannot represent. # - We **kept Glen A at the value of the first delivery** (calibrated there to the # ice volume measured from DEM differencing, 2021). This keeps all model parameters # identical, so that the change of divides is the only difference between the two # runs. The price is that the total volume is not forced to be the same. We could # have re-calibrated Glen A to match the same volume, but then the two runs would # differ by a parameter as well. # # The volume difference starts in 2025 (the inversion) and is carried back to 2000 # by the historical run. Relative changes (in % of 2020) are much less affected than # absolute volumes, which is why the comparisons below show both. # %% [markdown] # ## Helper functions # # Runoff is defined as in the first delivery (and the OGGM-EDU tutorial on # glacier water resources): on-glacier and off-glacier melt plus on-glacier and # off-glacier liquid precipitation, each clipped at zero, with an 11-year centred # smoothing. It is computed over a fixed area (the 2026 glacier outline), so it keeps # counting rain and snowmelt after the ice is gone. The hydro variables are only # stored in the projection files, so runoff starts in 2026. # %% ssps = ["ssp126", "ssp245", "ssp370", "ssp585"] ssp_colors = {"ssp126": "steelblue", "ssp245": "gold", "ssp370": "darkorange", "ssp585": "firebrick"} runoff_vars = ["melt_on_glacier", "melt_off_glacier", "liq_prcp_on_glacier", "liq_prcp_off_glacier"] def get_ids(version, icecap=None): """All rgi_ids of an inventory, or those of one ice cap.""" if icecap is None: return hist[version].rgi_id.values return icecap_ids[version][icecap] def regional_series(version, icecap=None): """Historical and projected regional totals for one inventory.""" ids = get_ids(version, icecap) h = hist[version].sel(rgi_id=ids) p = proj[version].sel(rgi_id=ids) mass_2020 = float(h.mass_kg.sum("rgi_id").sel(time=2020)) out = { "hist_volume": h.volume.sum("rgi_id") * 1e-9, # km³ "hist_area": h.area.sum("rgi_id") * 1e-6, # km² "hist_mass_pct": h.mass_kg.sum("rgi_id") / mass_2020 * 100, "volume": p.volume.sum("rgi_id") * 1e-9, "area": p.area.sum("rgi_id") * 1e-6, "mass_pct": p.mass_kg.sum("rgi_id") / mass_2020 * 100, } runoff = sum(p[v].clip(min=0) for v in runoff_vars).sum("rgi_id") * 1e-12 # Gt yr-1 out["runoff"] = runoff.rolling(time=11, center=True, min_periods=1).mean() return out def ssp_runs(da, ssp): return da.sel(gcm_scenario=[g for g in da.gcm_scenario.values if g.endswith(f"_{ssp}")]) def plot_comparison(icecap=None, fname=None): """4x4 figure: rows = volume, area, mass (% of 2020), runoff; columns = SSPs.""" series = {v: regional_series(v, icecap) for v in ["ice", "water"]} rows = [("volume", "Volume (km³)"), ("area", "Area (km²)"), ("mass_pct", "Mass (% of 2020)"), ("runoff", "Runoff (Gt yr$^{-1}$)")] fig, axes = plt.subplots(4, 4, figsize=(17, 14), sharex=True, sharey="row") for j, ssp in enumerate(ssps): c = ssp_colors[ssp] for i, (var, ylabel) in enumerate(rows): ax = axes[i, j] for version, ls in [("ice", "-"), ("water", "--")]: s = series[version] if f"hist_{var}" in s: hs = s[f"hist_{var}"].sel(time=slice(2000, 2026)) ax.plot(hs.time, hs, color="black", ls=ls, lw=1.2) da = ssp_runs(s[var], ssp) med = da.median("gcm_scenario") p10 = da.quantile(0.1, "gcm_scenario") p90 = da.quantile(0.9, "gcm_scenario") ax.plot(da.time, med, color=c, ls=ls, lw=1.8, label=f"{labels[version]}" if (i == 0 and j == 0) else None) if version == "ice": ax.fill_between(da.time, p10, p90, color=c, alpha=0.2, lw=0) else: ax.plot(da.time, p10, color=c, ls=":", lw=0.9) ax.plot(da.time, p90, color=c, ls=":", lw=0.9) if var == "mass_pct": ax.axhline(100, color="gray", lw=0.7, ls=":") ax.grid(True, alpha=0.3) if i == 0: n = len(ssp_runs(series["ice"]["volume"], ssp).gcm_scenario) ax.set_title(f"{ssp} (n={n})") if j == 0: ax.set_ylabel(ylabel) if i == 3: ax.set_xlabel("Year") axes[0, 0].legend(loc="lower left", fontsize=9) region = "Iceland" if icecap is None else icecap fig.suptitle(f"{region}: ice divides (solid, shaded 10-90%) vs water divides " f"(dashed, dotted 10-90%). Black: historical (CARRA).", fontsize=13) plt.tight_layout() if fname is not None: fig.savefig(f"{fig_dir}/{fname}", dpi=150, bbox_inches="tight") plt.show() return series def summary_table(series): """Median across GCMs of key numbers, per SSP and inventory.""" rows = [] for ssp in ssps: for version in ["ice", "water"]: s = series[version] vol = ssp_runs(s["volume"], ssp) mass = ssp_runs(s["mass_pct"], ssp) ro = ssp_runs(s["runoff"], ssp) rows.append({ "SSP": ssp, "Inventory": labels[version], "Volume 2020 (km³)": float(s["hist_volume"].sel(time=2020)), "Volume 2100 (km³)": float(vol.sel(time=2100).median()), "Mass 2050 (% of 2020)": float(mass.sel(time=2050).median()), "Mass 2100 (% of 2020)": float(mass.sel(time=2100).median()), "Peak water year": float(ro.idxmax("time").median()), "Peak runoff (Gt/yr)": float(ro.max("time").median()), }) return pd.DataFrame(rows).set_index(["SSP", "Inventory"]).round(1) # %% [markdown] # ## All of Iceland # %% series_iceland = plot_comparison(fname="comparison_iceland.png") summary_table(series_iceland) # %% [markdown] # For Iceland as a whole, the two inventories are very close: # # - The water-divide run has about 2% less ice in 2020 (3220 vs 3276 km³): the ice # thickness inversion depends on the shape of the entities. # - Relative to 2020, it loses nearly the same fraction of its mass by 2100, slightly # more under high emissions (32.8% vs 34.7% left under ssp585). # - Runoff is 2-3% higher at its peak. The peak-water year is unchanged except under # ssp126 (2051 vs 2045), where the runoff curve is so flat that the year of the # maximum is not robust. # # All these differences are much smaller than the spread across climate models. # # The differences come from the ice caps that were re-divided. All other glaciers # give identical results in both runs. # %% [markdown] # ## Ice caps most relevant for hydropower # # We repeat the comparison for three ice caps: # # - **Vatnajökull**: by far the largest ice mass; feeds Hálslón (Kárahnjúkar) and the # Tungnaá branch of the Þjórsá system. # - **Hofsjökull**: headwaters of Þjórsá and Blanda. # - **Langjökull**: where the change of divides has the largest effect; the # Blanda catchment entity of the inventory lies on Langjökull. # # The same function works for any ice cap in `data/inventory/icecap_ids_rgi2025hydro.json` # (Mýrdalsjökull, Drangajökull, ...). # %% series_vatna = plot_comparison("Vatnajökull", fname="comparison_vatnajokull.png") summary_table(series_vatna) # %% series_hofs = plot_comparison("Hofsjökull", fname="comparison_hofsjokull.png") summary_table(series_hofs) # %% series_lang = plot_comparison("Langjökull", fname="comparison_langjokull.png") summary_table(series_lang) # %% [markdown] # At the ice-cap level the differences are larger than for Iceland as a whole, but # remain moderate compared to the climate-model spread: # # - **Vatnajökull**: about 3% less ice in the water-divide run (2717 vs 2798 km³ in # 2020), but nearly identical relative mass loss, runoff and peak-water timing. # - **Hofsjökull**: about 1% more ice, a slightly faster relative mass loss # (about 3 percentage points less left in 2100), and a runoff peak that is # slightly higher and equal or a few years later. # - **Langjökull**: the largest change. The water-divide run has about 14% more ice # (188 vs 165 km³; fewer, larger entities give thicker ice in the inversion). # It keeps slightly more of its 2020 mass until about 2050, then loses it faster # (e.g. 13% vs 18% left in 2100 under ssp245). Its runoff peak is 10-13% higher and # 3-10 years later in all scenarios. # # Note that the total glacier area of each ice cap is identical in both inventories, # so the differences come from how the ice is distributed between entities. For # individual catchments, see notebook `02_catchments`.