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/.

In [1]:
import json

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

Load data¶

In [2]:
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)"}

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.

In [3]:
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")

69 CMIP6 runs used in both inventories

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.

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.

In [4]:
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)

All of Iceland¶

In [5]:
series_iceland = plot_comparison(fname="comparison_iceland.png")
summary_table(series_iceland)
No description has been provided for this image
Out[5]:
Volume 2020 (km³) Volume 2100 (km³) Mass 2050 (% of 2020) Mass 2100 (% of 2020) Peak water year Peak runoff (Gt/yr)
SSP Inventory
ssp126 Ice divides (RGI2025) 3275.6 1894.6 83.8 57.9 2045.0 41.4
Water divides (hydro) 3219.9 1858.8 83.8 57.8 2051.0 42.4
ssp245 Ice divides (RGI2025) 3275.6 1759.0 83.0 53.8 2055.0 44.9
Water divides (hydro) 3219.9 1703.6 83.0 53.0 2055.0 45.8
ssp370 Ice divides (RGI2025) 3275.6 1454.9 83.5 44.5 2072.0 47.1
Water divides (hydro) 3219.9 1387.6 83.5 43.2 2072.0 48.3
ssp585 Ice divides (RGI2025) 3275.6 1133.4 79.5 34.7 2072.0 52.7
Water divides (hydro) 3219.9 1054.3 79.4 32.8 2072.0 53.7

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.

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, ...).

In [6]:
series_vatna = plot_comparison("Vatnajökull", fname="comparison_vatnajokull.png")
summary_table(series_vatna)
No description has been provided for this image
Out[6]:
Volume 2020 (km³) Volume 2100 (km³) Mass 2050 (% of 2020) Mass 2100 (% of 2020) Peak water year Peak runoff (Gt/yr)
SSP Inventory
ssp126 Ice divides (RGI2025) 2798.0 1726.0 86.3 61.7 2052.0 31.2
Water divides (hydro) 2716.8 1676.5 86.3 61.8 2052.0 31.2
ssp245 Ice divides (RGI2025) 2798.0 1637.0 85.3 58.6 2060.5 32.4
Water divides (hydro) 2716.8 1591.0 85.3 58.7 2058.0 32.9
ssp370 Ice divides (RGI2025) 2798.0 1364.1 86.1 48.8 2072.0 35.5
Water divides (hydro) 2716.8 1291.4 86.2 47.7 2072.0 36.6
ssp585 Ice divides (RGI2025) 2798.0 1080.1 82.6 38.7 2076.0 43.1
Water divides (hydro) 2716.8 1018.0 82.5 37.6 2076.0 42.7
In [7]:
series_hofs = plot_comparison("Hofsjökull", fname="comparison_hofsjokull.png")
summary_table(series_hofs)
No description has been provided for this image
Out[7]:
Volume 2020 (km³) Volume 2100 (km³) Mass 2050 (% of 2020) Mass 2100 (% of 2020) Peak water year Peak runoff (Gt/yr)
SSP Inventory
ssp126 Ice divides (RGI2025) 160.8 59.5 73.5 36.9 2045.0 3.1
Water divides (hydro) 162.5 55.9 73.2 34.3 2052.0 3.3
ssp245 Ice divides (RGI2025) 160.8 51.5 74.3 32.0 2050.0 3.2
Water divides (hydro) 162.5 47.3 73.9 29.1 2052.5 3.4
ssp370 Ice divides (RGI2025) 160.8 37.6 74.1 23.3 2053.0 3.5
Water divides (hydro) 162.5 32.4 73.8 19.9 2057.0 3.7
ssp585 Ice divides (RGI2025) 160.8 20.9 69.0 13.0 2053.0 3.7
Water divides (hydro) 162.5 16.3 68.2 10.0 2053.0 4.2
In [8]:
series_lang = plot_comparison("Langjökull", fname="comparison_langjokull.png")
summary_table(series_lang)
No description has been provided for this image
Out[8]:
Volume 2020 (km³) Volume 2100 (km³) Mass 2050 (% of 2020) Mass 2100 (% of 2020) Peak water year Peak runoff (Gt/yr)
SSP Inventory
ssp126 Ice divides (RGI2025) 164.7 36.0 66.8 21.6 2044.0 3.8
Water divides (hydro) 187.5 31.3 68.5 16.5 2054.0 4.3
ssp245 Ice divides (RGI2025) 164.7 30.3 67.4 18.3 2050.0 3.9
Water divides (hydro) 187.5 25.2 69.0 13.3 2053.5 4.3
ssp370 Ice divides (RGI2025) 164.7 15.7 69.8 9.5 2051.0 4.2
Water divides (hydro) 187.5 11.8 71.4 6.2 2056.0 4.7
ssp585 Ice divides (RGI2025) 164.7 3.5 61.8 2.1 2051.0 4.7
Water divides (hydro) 187.5 2.1 63.6 1.1 2055.0 5.3

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.