Using the water-divide runs: catchments

This notebook shows how to work with the new inventory (RGI2025-hydro v2) at the level of individual hydrological catchments, which is what matters for hydropower. It gives:

1. the structure of the inventory (which entity is what);
2. a summary table for all catchments (volume, mass loss, peak water);
3. a worked example for one aggregate catchment, Hálslón.

All projections are CMIP6 to 2100 (median and 10th-90th percentile range across GCMs). The same runs also go to 2300 and exist for CMIP5: see data/projections/.

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

data_dir = "data"
fig_dir = "figures"
table_dir = "tables"

The inventory

The inventory data/inventory/RGI2025hydro-v7.0-G-06_iceland.shp follows the RGI7 format, with a few extra columns:

• ent_type: catchment (one of your water-divide catchments), remnant (ice not covered by any catchment, modelled so that totals stay consistent, but not assigned to a river), or rgi2025 (unchanged glacier from the first delivery, with the same rgi_id).
• ws_id and catchment: the ID and name of your water-divide polygon.
• icecap: the ice cap the entity belongs to.

New entities have IDs from RGI2025-v7.0-G-06-10001 onwards, so they never collide with RGI2025 IDs.

inv = gpd.read_file(f"{data_dir}/inventory/RGI2025hydro-v7.0-G-06_iceland.shp")
inv = inv.set_index("rgi_id")
inv.groupby("ent_type").area_km2.agg(["count", "sum"]).round(1)

ERROR 1: PROJ: proj_create_from_database: Open of /home/users/fmaussion/.miniforge3/envs/oggm_env/share/proj failed

	count	sum
ent_type
catchment	59	9547.6
remnant	20	219.9
rgi2025	492	431.8

Maps of the new entities, ice cap by ice cap, are in figures/inventory_maps/. Your original water divides are drawn as a red dashed line on top of the modelled entities.

From your catchments to model entities

Most catchments correspond to exactly one model entity (the one with the same ws_id). Two of your catchments contain smaller ones and are modelled as sums (data/inventory/catchment_lookup.csv):

• Halslon10 (0035) = Kring + illakv + jokulsa + jokulkvisl;
• Jökulsá í Lóni (0041) = Skyndidalsá + the rest of the Jökulsá í Lóni catchment.

Þingvallavatn is modelled with the Þingvallavatn_2 delineation (0023), as you chose. It includes 43 km² that are also part of your Hagavatn polygon: in the model, this area belongs to Þingvallavatn only.

lookup = pd.read_csv(f"{data_dir}/inventory/catchment_lookup.csv", dtype={"ws_id": str})

catchment_ids = {}  # catchment name -> list of rgi_ids
for _, row in lookup.iterrows():
    if row.ws_id in ["0035", "0041"]:  # aggregates: sum of sub-catchments
        catchment_ids[row.catchment] = row.rgi_ids.split(";")
    elif row.ws_id in inv.ws_id.values:  # one entity per catchment
        catchment_ids[row.catchment] = list(inv.index[inv.ws_id == row.ws_id])
print(f"{len(catchment_ids)} catchments (incl. 2 aggregates)")
print("Hálslón:", catchment_ids["Halslon10"])

60 catchments (incl. 2 aggregates)
Hálslón: ['RGI2025-v7.0-G-06-10041', 'RGI2025-v7.0-G-06-10047', 'RGI2025-v7.0-G-06-10048', 'RGI2025-v7.0-G-06-10049']

Load the runs

hist = xr.open_dataset(f"{data_dir}/historical/hydro_v2/run_output_spinup_historical.nc")
proj = xr.open_dataset(f"{data_dir}/projections/hydro_v2/merged_2100_projections_cmip6.nc")

# See notebook 01: four runs lose many small glaciers and are excluded here
exclude = ["CESM2-WACCM_ssp126", "CESM2-WACCM_ssp534-over",
           "MRI-ESM2-0_ssp126", "MRI-ESM2-0_ssp534-over"]
proj = proj.sel(gcm_scenario=[g for g in proj.gcm_scenario.values if g not in exclude])

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 ssp_runs(da, ssp):
    return da.sel(gcm_scenario=[g for g in da.gcm_scenario.values if g.endswith(f"_{ssp}")])


def catchment_series(ids):
    """Mass (% of 2020) and smoothed annual runoff (Gt/yr) for a list of entities."""
    h = hist.sel(rgi_id=ids)
    p = proj.sel(rgi_id=ids)
    mass_2020 = float(h.mass_kg.sum("rgi_id").sel(time=2020))
    runoff = sum(p[v].clip(min=0) for v in runoff_vars).sum("rgi_id") * 1e-12
    return {
        "volume_2020": float(h.volume.sum("rgi_id").sel(time=2020)) * 1e-9,
        "hist_mass_pct": h.mass_kg.sum("rgi_id") / mass_2020 * 100,
        "mass_pct": p.mass_kg.sum("rgi_id") / mass_2020 * 100,
        "runoff": runoff.rolling(time=11, center=True, min_periods=1).mean(),
    }

Summary table for all catchments

Median across GCMs. "Peak water" is the year of maximum 11-yr smoothed runoff. For small catchments the GCM spread is large: use the full time series (next section) rather than single numbers when possible.

rows = []
for name, ids in catchment_ids.items():
    s = catchment_series(ids)
    row = {
        "Catchment": name,
        "Ice cap": inv.loc[ids[0], "icecap"],
        "Area 2025 (km²)": inv.loc[ids, "area_km2"].sum(),
        "Volume 2020 (km³)": s["volume_2020"],
    }
    for ssp in ["ssp245", "ssp585"]:
        mass = ssp_runs(s["mass_pct"], ssp)
        ro = ssp_runs(s["runoff"], ssp)
        row[f"Mass 2050 {ssp} (%)"] = float(mass.sel(time=2050).median())
        row[f"Mass 2100 {ssp} (%)"] = float(mass.sel(time=2100).median())
        row[f"Peak water {ssp}"] = float(ro.idxmax("time").median())
    rows.append(row)

summary = pd.DataFrame(rows).set_index("Catchment").sort_values(["Ice cap", "Area 2025 (km²)"])
summary.to_csv(f"{table_dir}/catchment_summary.csv")
summary.round(1)

	Ice cap	Area 2025 (km²)	Volume 2020 (km³)	Mass 2050 ssp245 (%)	Mass 2100 ssp245 (%)	Peak water ssp245	Mass 2050 ssp585 (%)	Mass 2100 ssp585 (%)	Peak water ssp585
Catchment
Tharalatursos	Drangajökull	6.6	0.5	32.4	0.0	2029.5	18.7	0.0	2039.0
Bjarnafjardara	Drangajökull	16.2	1.4	39.3	0.4	2041.0	25.3	0.0	2042.0
Reykjafjardaros	Drangajökull	18.0	1.8	46.6	0.8	2038.0	33.9	0.0	2048.0
Sela	Drangajökull	25.4	3.0	70.4	6.6	2058.0	63.2	0.0	2056.0
Jokulsa_i_leirufirdi	Drangajökull	26.4	3.0	65.7	1.2	2050.0	57.0	0.0	2056.0
Morilla	Drangajökull	33.9	4.5	65.3	5.9	2050.0	58.2	0.0	2056.0
Hofsjökull-vatnasvið-2008-sv_e	Hofsjökull	0.7	0.1	35.8	0.0	2040.0	28.8	0.0	2042.0
Hofsjökull-vatnasvið-2008-sv_b	Hofsjökull	2.0	0.3	46.8	0.5	2040.5	40.6	0.0	2044.0
Hofsjökull-vatnasvið-2008-sv_c	Hofsjökull	4.1	0.7	65.6	3.6	2053.5	59.5	0.0	2054.0
Hofsjökull-vatnasvið-2008-sv_d	Hofsjökull	5.9	0.9	59.2	0.3	2050.0	53.3	0.0	2052.0
Hofsjökull-vatnasvið-2008-sv_g	Hofsjökull	25.0	5.2	74.7	19.8	2065.0	69.9	0.6	2069.0
Hofsjökull-vatnasvið-2008-sv_r	Hofsjökull	45.4	9.0	68.5	17.6	2052.5	64.1	4.7	2055.0
Hofsjökull-vatnasvið-2008-sv_f	Hofsjökull	50.3	12.8	85.5	30.5	2079.0	79.0	4.9	2071.0
Hofsjökull-vatnasvið-2008-sv_q	Hofsjökull	56.6	13.8	74.4	26.0	2061.5	69.4	6.7	2060.0
Hofsjökull-vatnasvið-2008-sv_s	Hofsjökull	106.1	24.3	78.3	36.7	2057.0	73.5	18.6	2055.0
Hofsjökull-vatnasvið-2008-sv_a	Hofsjökull	164.0	32.5	81.0	42.5	2058.5	75.2	21.5	2053.0
Hofsjökull-vatnasvið-2008-sv_p	Hofsjökull	179.3	38.3	62.5	19.0	2049.0	55.6	4.2	2051.0
Blanda	Langjökull	11.8	1.3	41.2	0.0	2045.0	32.0	0.0	2045.0
Hagavatn	Langjökull	125.1	33.6	73.8	29.1	2052.5	69.0	5.9	2052.0
Þingvallavatn_2	Langjökull	125.5	29.9	71.7	16.8	2055.5	66.0	0.0	2055.0
Hvíta_Borgarfirði	Langjökull	275.1	53.5	63.2	3.0	2053.0	56.9	0.0	2053.0
Hvítá	Langjökull	276.7	69.3	70.6	12.5	2057.5	65.8	0.2	2062.0
Krossá	Mýrdalsjökull	21.9	3.1	70.8	37.1	2029.5	66.1	12.1	2036.0
Klifandi_Mýrdalur	Mýrdalsjökull	35.9	3.9	51.2	22.0	2026.5	45.1	2.7	2027.0
Bláfjallakvísl_Bláfellsá_Jökulkvísl	Mýrdalsjökull	36.0	9.8	67.7	29.8	2033.5	62.9	14.6	2046.0
Hólmsá	Mýrdalsjökull	46.1	10.0	49.4	4.6	2029.5	46.6	0.1	2043.0
Leirá	Mýrdalsjökull	53.4	8.9	55.0	33.4	2028.5	49.3	6.9	2032.0
Skóga_og_Sólheimasandar	Mýrdalsjökull	64.2	14.9	73.9	44.3	2031.5	71.0	16.5	2064.0
Markarfljót_án_Krossár	Mýrdalsjökull	113.7	26.8	69.0	33.8	2031.0	66.5	17.8	2043.0
Mýrdalssandur	Mýrdalsjökull	138.2	41.4	78.4	55.7	2033.0	75.0	36.7	2050.0
v_skaftarketill	Vatnajökull	25.0	6.2	86.2	33.3	2089.0	81.5	0.0	2079.0
Brunná	Vatnajökull	28.6	10.3	55.0	0.0	2056.5	48.8	0.0	2052.0
kverka	Vatnajökull	32.4	8.7	36.2	0.0	2036.0	32.4	0.0	2038.0
skyndidalsá_	Vatnajökull	33.6	5.6	53.5	15.2	2036.0	48.8	1.4	2043.0
sylgja	Vatnajökull	37.4	7.7	57.1	0.1	2045.5	48.9	0.0	2051.0
e_skaftarketill	Vatnajökull	39.3	8.2	88.1	42.5	2087.0	83.5	0.4	2077.0
jokulkvisl	Vatnajökull	40.1	11.1	79.4	20.3	2068.5	75.3	0.0	2071.0
Djúpá	Vatnajökull	69.8	18.0	64.3	3.0	2053.0	57.7	0.0	2052.0
Jökulsá_í_lóni_	Vatnajökull	93.2	14.5	55.1	10.5	2040.5	50.3	0.9	2044.0
Öræfajökull_vatnaskila_butar	Vatnajökull	93.9	22.7	79.6	51.3	2059.0	75.3	9.9	2074.0
skjalfandfljot_	Vatnajökull	98.8	14.9	84.1	37.5	2067.0	79.7	12.2	2070.0
tungnaa	Vatnajökull	101.5	31.8	67.6	11.4	2047.5	61.6	0.0	2051.0
Jokfljotsdvsk	Vatnajökull	121.6	25.6	79.9	24.5	2057.5	76.3	0.8	2067.0
hólmsá	Vatnajökull	156.9	45.9	78.1	45.4	2046.5	73.5	22.9	2052.0
grimsvotn_	Vatnajökull	173.2	37.2	106.9	111.6	2077.5	104.3	49.4	2086.0
Heinabergsvötn_	Vatnajökull	221.2	66.5	77.2	39.1	2053.5	72.9	16.6	2055.0
Breiðá_Fjallsá_	Vatnajökull	225.4	58.4	69.1	53.2	2028.0	67.4	33.7	2032.0
Hornafjardarfljót_	Vatnajökull	228.4	67.8	77.3	38.8	2047.5	72.2	5.4	2074.0
Súla_	Vatnajökull	248.7	92.3	76.9	26.4	2067.0	73.8	9.5	2062.0
Kreppa	Vatnajökull	286.6	87.1	91.6	69.5	2092.5	89.3	46.6	2090.0
Hverfisfljot_	Vatnajökull	303.6	131.5	88.3	54.1	2086.5	85.2	24.5	2087.0
jokulsa	Vatnajökull	339.1	111.1	88.0	50.8	2083.0	84.9	19.2	2083.0
kaldakvisl	Vatnajökull	347.4	105.1	85.3	48.0	2068.5	81.8	26.1	2076.0
illakv	Vatnajökull	358.1	122.7	93.2	66.9	2090.5	90.8	39.1	2090.0
skafta	Vatnajökull	375.3	154.6	81.8	51.8	2054.0	78.3	31.3	2076.0
Kring	Vatnajökull	443.6	163.1	94.3	77.4	2084.5	92.1	55.3	2093.0
Jokulsa_a_Breidam_sandi_	Vatnajökull	719.1	292.4	81.5	56.5	2049.0	78.7	38.6	2063.0
jokulsa_a_fjollum_	Vatnajökull	1115.9	452.4	92.2	69.8	2090.5	90.2	44.4	2087.0
Skeiðará_	Vatnajökull	1133.8	521.4	87.1	62.2	2067.5	84.5	44.1	2071.0
Halslon10	Vatnajökull	1181.0	408.0	91.9	65.8	2091.0	89.3	39.1	2083.0

The table is also saved as tables/catchment_summary.csv.

Caveat: catchments in the accumulation area

OGGM treats each catchment as a glacier with its own flowline. In reality, ice in the interior of an ice cap flows across water divides towards the outlet glaciers. A catchment that lies mostly in the accumulation area therefore has no realistic way to lose its ice in the model, and may even thicken. These catchments stand out in the table: they keep almost all of their mass until 2050, and some gain mass.

interior = summary[summary["Mass 2050 ssp245 (%)"] > 90]
interior[["Ice cap", "Area 2025 (km²)", "Mass 2050 ssp245 (%)",
          "Mass 2100 ssp245 (%)", "Peak water ssp245"]].round(1)

	Ice cap	Area 2025 (km²)	Mass 2050 ssp245 (%)	Mass 2100 ssp245 (%)	Peak water ssp245
Catchment
grimsvotn_	Vatnajökull	173.2	106.9	111.6	2077.5
Kreppa	Vatnajökull	286.6	91.6	69.5	2092.5
illakv	Vatnajökull	358.1	93.2	66.9	2090.5
Kring	Vatnajökull	443.6	94.3	77.4	2084.5
jokulsa_a_fjollum_	Vatnajökull	1115.9	92.2	69.8	2090.5
Halslon10	Vatnajökull	1181.0	91.9	65.8	2091.0

grimsvotn is the clearest case: it gains mass throughout the century under ssp245, which is not realistic. The other catchments in this list are the large northern Vatnajökull catchments, including all four that make up Hálslón (below). They extend high into the accumulation area, so a slower mass loss and a late peak water are partly expected, but part of it is likely this artefact, and the runs alone cannot separate the two.

For these catchments, treat the mass and peak-water numbers with caution. Their runoff in the coming decades is more robust, because it is driven by the surface melt and rain in the catchment rather than by ice flow.

Example: Hálslón

The Hálslón catchment (Halslon10) is the sum of four model entities. Any other catchment, or any combination of entities, works the same way: replace ids.

ids = catchment_ids["Halslon10"]
s = catchment_series(ids)
print(f"Area {inv.loc[ids, 'area_km2'].sum():.0f} km², volume 2020 {s['volume_2020']:.0f} km³")

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
hs = s["hist_mass_pct"].sel(time=slice(2000, 2026))
axes[0].plot(hs.time, hs, color="black", lw=1.5, label="Historical (CARRA)")
for ssp in ssps:
    for ax, var in zip(axes, ["mass_pct", "runoff"]):
        da = ssp_runs(s[var], ssp)
        c = ssp_colors[ssp]
        ax.plot(da.time, da.median("gcm_scenario"), color=c,
                label=f"{ssp} (n={da.gcm_scenario.size})")
        ax.fill_between(da.time, da.quantile(0.1, "gcm_scenario"),
                        da.quantile(0.9, "gcm_scenario"), color=c, alpha=0.2, lw=0)
axes[0].axhline(100, color="gray", lw=0.7, ls=":")
axes[0].set_ylabel("Mass (% of 2020)")
axes[1].set_ylabel("Annual runoff (Gt yr$^{-1}$, 11-yr smoothing)")
for ax in axes:
    ax.set_xlabel("Year")
    ax.grid(True, alpha=0.3)
    ax.legend(fontsize=9)
fig.suptitle("Hálslón catchment (Kring + illakv + jokulsa + jokulkvisl)", fontsize=13)
plt.tight_layout()
fig.savefig(f"{fig_dir}/halslon_mass_runoff.png", dpi=150, bbox_inches="tight")
plt.show()

Area 1181 km², volume 2020 408 km³

Seasonal cycle of runoff

The projection files also store monthly runoff (variables ending in _monthly, in kg per month). Here: the mean seasonal cycle of Hálslón runoff for three decades, under ssp245 (median across GCMs). Runoff is converted to an average discharge in m³ s⁻¹.

ssp = "ssp245"
p = ssp_runs(proj.sel(rgi_id=ids), ssp)
monthly = sum(p[f"{v}_monthly"].clip(min=0) for v in runoff_vars).sum("rgi_id")
seconds_per_month = 365 * 24 * 3600 / 12  # the model uses a 365-day calendar
discharge = monthly / 1000 / seconds_per_month  # kg/month -> m³/s

decades = {"2026-2035": "lightsteelblue", "2056-2065": "steelblue", "2091-2100": "navy"}
fig, ax = plt.subplots(figsize=(9, 5))
for (period, c) in decades.items():
    y0, y1 = [int(y) for y in period.split("-")]
    cycle = discharge.sel(time=slice(y0, y1)).mean("time").median("gcm_scenario")
    ax.plot(p.calendar_month_2d, cycle, color=c, marker="o", label=period)
ax.set_xticks(range(1, 13))
ax.set_xticklabels(["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"])
ax.set_ylabel("Mean discharge (m³ s$^{-1}$)")
ax.set_xlabel("Month")
ax.grid(True, alpha=0.3)
ax.legend()
ax.set_title(f"Hálslón catchment: seasonal cycle of glacier runoff ({ssp}, median)")
plt.tight_layout()
fig.savefig(f"{fig_dir}/halslon_seasonal_cycle.png", dpi=150, bbox_inches="tight")
plt.show()

Note: as in the first delivery, runoff is computed over a fixed area (the glacier outline at the start of the projections). It therefore includes rain and snowmelt on terrain that becomes ice-free, but not the rest of the river catchment.