{ "cells": [ { "cell_type": "markdown", "id": "9df71dc0", "metadata": {}, "source": [ "# Iceland OGGM — Experiment Selection & Validation\n", "\n", "Before running the CMIP5/CMIP6 projections, we calibrated OGGM against several\n", "independent observations (glacier area 2010/2025, Hugonnet et al. 2000-2020 mass\n", "change, DEM-differenced volume 2021, and Andri's mass-balance maps).\n", "This was done across a sweep of configuration choices:\n", "\n", "- **climate forcing**: ERA5 vs. CARRA reanalysis (1991-2026)\n", "- **flowline bed shape**: trapezoidal, with the shape parameter lambda ranging from\n", " 2 (~45° angle of side wall, measured from horizontal) to 16 (~7° angle, much flatter\n", " than 45°)\n", "- **mass-balance models**: monthly temperature-index model with surface type tracking\n", " (default) and without tracking (Monthly TI)\n", "\n", "**Outcome: CARRA, trapezoidal lambda = 14 was chosen as the delivered configuration.**\n", "For ERA5, lambda = 8 was the closest contender and is shown throughout for comparison —\n", "the difference between the two is small. The Monthly TI variant (without snow tracking)\n", "of each was also tested and is included alongside its parent configuration.\n", "\n", "**The choice is a compromise, not a perfect fit.** The lambda sweep (Appendix) shows a\n", "systematic trade-off: increasing lambda improves the area match (particularly the 2010\n", "outline, an independent check — 2025 is matched by construction, since it is the\n", "spinup's initialisation target) but pushes the modelled 2000-2020 volume loss further\n", "from the Hugonnet estimate; decreasing lambda does the opposite. No single lambda\n", "simultaneously minimises both mismatches. CARRA lambda=14 and ERA5 lambda=8 are the\n", "points in the sweep that balance the two best — not the configurations that match\n", "either observation alone most closely.\n", "\n", "This notebook first summarises that decision (see *Decision summary* below), then walks\n", "through the validation evidence for the shortlisted configurations. The full lambda\n", "sweep that led to the shortlist is kept in the **Appendix** at the end, for transparency\n", "and reproducibility — it is not part of the main narrative.\n", "\n", "**On data availability:** to keep the delivery package simple and clear, only the\n", "NetCDF output of the delivered configuration (and, for comparison, the ERA5\n", "alternative) ships alongside this notebook — not the full ~15-configuration sweep shown\n", "in the Appendix. That full sweep's data can be made available on request. As shipped,\n", "this notebook reads from `/home/shared/fmaussion_pschmitt/iceland_projections` on the\n", "OGGM cluster and will **not run as-is** without access to that path — it is included to\n", "document the calibration process and decision, not as an executable artefact for the\n", "recipient." ] }, { "cell_type": "code", "execution_count": 1, "id": "4137881f", "metadata": { "execution": { "iopub.execute_input": "2026-07-07T21:44:04.646417Z", "iopub.status.busy": "2026-07-07T21:44:04.646134Z", "iopub.status.idle": "2026-07-07T21:44:05.673621Z", "shell.execute_reply": "2026-07-07T21:44:05.672570Z" } }, "outputs": [], "source": [ "import os\n", "from IPython.display import HTML\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "import yaml\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "b5f0ba4b", "metadata": {}, "source": [ "## Load model output and observations\n", "\n", "`CONFIGS` lists every configuration in the calibration sweep (label -> working\n", "directory suffix). `SHORTLIST` is the subset shown in the main narrative below; the\n", "rest only reappears in the Appendix." ] }, { "cell_type": "code", "execution_count": 2, "id": "6749dd56", "metadata": { "execution": { "iopub.execute_input": "2026-07-07T21:44:05.676353Z", "iopub.status.busy": "2026-07-07T21:44:05.676070Z", "iopub.status.idle": "2026-07-07T21:44:06.734451Z", "shell.execute_reply": "2026-07-07T21:44:06.733189Z" } }, "outputs": [], "source": [ "project_dir = \"/home/shared/fmaussion_pschmitt/iceland_projections\"\n", "workflow_dir = f\"{project_dir}/complete_workflow\"\n", "validation_data_dir = f\"{project_dir}/validation/validation_data\"\n", "\n", "fig_dir = \"figures\" # relative to packaging_delivery/, kept self-contained for delivery\n", "os.makedirs(fig_dir, exist_ok=True)\n", "\n", "CONFIGS = {\n", " \"era5_2\": \"ERA5, trapez. λ=2 (~45°)\",\n", " \"era5_4\": \"ERA5, trapez. λ=4 (~26°)\",\n", " \"era5_8\": \"ERA5, trapez. λ=8 (~14°)\",\n", " \"era5_14\": \"ERA5, trapez. λ=14 (~8°)\",\n", " \"era5_16\": \"ERA5, trapez. λ=16 (~7°)\",\n", " \"era5_8_MonthlyTI\": \"ERA5, trapez. λ=8 (~14°), Monthly TI\",\n", " \"era5_14_MonthlyTI\": \"ERA5, trapez. λ=14 (~8°), Monthly TI\",\n", " \"carra_2\": \"CARRA, trapez. λ=2 (~45°)\",\n", " \"carra_4\": \"CARRA, trapez. λ=4 (~26°)\",\n", " \"carra_8\": \"CARRA, trapez. λ=8 (~14°)\",\n", " \"carra_12\": \"CARRA, trapez. λ=12 (~9°)\",\n", " \"carra_14\": \"CARRA, trapez. λ=14 (~8°)\",\n", " \"carra_16\": \"CARRA, trapez. λ=16 (~7°)\",\n", " \"carra_8_MonthlyTI\": \"CARRA, trapez. λ=8 (~14°), Monthly TI\",\n", " \"carra_14_MonthlyTI\": \"CARRA, trapez. λ=14 (~8°), Monthly TI\",\n", "}\n", "\n", "# Delivered configuration (carra_14) + closest contender (era5_8), each with its\n", "# Monthly TI variant. Everything else in CONFIGS is sweep-only, shown in the Appendix.\n", "SHORTLIST = [\"carra_14\", \"carra_14_MonthlyTI\", \"era5_8\", \"era5_8_MonthlyTI\"]\n", "\n", "models = {\n", " key: xr.open_dataset(\n", " os.path.join(workflow_dir, f\"working_dir_adapted_{key}\",\n", " \"run_output_spinup_historical.nc\"))\n", " for key in CONFIGS\n", "}\n", "\n", "# area and volume-change observations, provided by Andri\n", "with open(os.path.join(validation_data_dir, \"area_dV_obs.yaml\")) as f:\n", " area_dv_obs = yaml.safe_load(f)\n", "\n", "# specific mass-balance observations derived from Andri's SMB maps\n", "ds_smb = xr.open_dataset(os.path.join(validation_data_dir, \"smbs_from_maps.nc\"))\n", "ds_smb[\"time_int\"] = ds_smb.time_a.dt.year" ] }, { "cell_type": "markdown", "id": "a72da947", "metadata": {}, "source": [ "## Helper functions\n", "\n", "These operate on the `models` dict and are reused for both the shortlist plots below\n", "and the full-sweep plots in the Appendix, so the two only differ in which keys they\n", "are called with." ] }, { "cell_type": "code", "execution_count": 3, "id": "a141bfd9", "metadata": { "execution": { "iopub.execute_input": "2026-07-07T21:44:06.738068Z", "iopub.status.busy": "2026-07-07T21:44:06.737737Z", "iopub.status.idle": "2026-07-07T21:44:06.777675Z", "shell.execute_reply": "2026-07-07T21:44:06.776697Z" }, "lines_to_next_cell": 1 }, "outputs": [], "source": [ "rgi_ids_volume_ddem = area_dv_obs[\"rgi_ids_volume_change_2010_2021\"]\n", "\n", "\n", "def plot_area(ax, keys):\n", " for key in keys:\n", " (models[key].sum(dim=\"rgi_id\").area / 1e6).plot(ax=ax, label=CONFIGS[key])\n", " ax.plot(2010, area_dv_obs[\"area_2010_km2\"], \"o\", color=\"k\", label=\"Outline area ~2010\")\n", " ax.plot(2025, area_dv_obs[\"area_2025_km2\"], \"o\", color=\"0.4\", label=\"Outline area ~2025\")\n", " ax.set_ylabel(\"Area (km²)\")\n", " ax.grid(True)\n", " ax.legend(bbox_to_anchor=[1.01, 0.5], loc=\"center left\")\n", "\n", "\n", "def plot_mass_change(ax, keys):\n", " for key in keys:\n", " ds = models[key].sum(dim=\"rgi_id\")\n", " (ds.mass_kg - ds.mass_kg.sel(time=2000)).plot(ax=ax, label=CONFIGS[key])\n", " var, err_var = \"hugonnet_adapted_2000_2020_kg\", \"hugonnet_2000-01-01_2020-01-01_kg_err\"\n", " ax.plot([2000, 2020], [0, area_dv_obs[var]], c=\"k\")\n", " ax.plot([2020, 2020],\n", " [area_dv_obs[var] - area_dv_obs[err_var], area_dv_obs[var] + area_dv_obs[err_var]],\n", " c=\"k\", label=\"Hugonnet (adapted) 2000-2020 ± error\")\n", " ax.set_ylabel(\"Mass change since 2000 (kg)\")\n", " ax.grid(True)\n", " ax.legend(bbox_to_anchor=[1.01, 0.5], loc=\"center left\")\n", "\n", "\n", "def plot_volume_bed(ax, keys):\n", " for key in keys:\n", " models[key].sel(rgi_id=rgi_ids_volume_ddem).volume.sum(dim=\"rgi_id\").plot(\n", " ax=ax, label=CONFIGS[key])\n", " ax.plot(2010, area_dv_obs[\"volume_2010\"], \"o\", color=\"k\", label=\"Volume 2010 (DEM - bed)\")\n", " ax.plot(2021, area_dv_obs[\"volume_2021\"], \"o\", color=\"0.4\", label=\"Volume 2021 (DEM - bed)\")\n", " ax.set_ylabel(\"Volume (m³)\")\n", " ax.grid(True)\n", " ax.legend(bbox_to_anchor=[1.01, 0.5], loc=\"center left\")\n", "\n", "\n", "def plot_volume_change_ddem(ax, keys):\n", " for key in keys:\n", " ds_vol = models[key].sel(rgi_id=rgi_ids_volume_ddem).volume.sum(dim=\"rgi_id\")\n", " (ds_vol - ds_vol.sel(time=2010)).plot(ax=ax, label=CONFIGS[key])\n", " ax.plot([2010, 2021], [0, area_dv_obs[\"volume_change_2010_2021\"]], c=\"k\",\n", " label=\"Volume change from DEM differencing\")\n", " ax.set_ylabel(\"Volume change since 2010 (m³)\")\n", " ax.grid(True)\n", " ax.legend(bbox_to_anchor=[1.01, 0.5], loc=\"center left\")\n", "\n", "\n", "def regional_smb(ds):\n", " \"\"\"Area-summed specific mass balance (mm w.e. yr-1) for the glaciers with SMB maps.\"\"\"\n", " ds = ds.sel(rgi_id=ds_smb.rgi_id).sum(dim=\"rgi_id\")\n", " return ds.volume.diff(dim=\"time\", label=\"lower\").reindex(time=ds.time) / ds.area * 900\n", "\n", "\n", "# Area weights for the regionally-aggregated observed SMB. Using the delivered\n", "# CARRA λ=14 config's glacier areas as the reference; areas barely differ between\n", "# configs at this stage, so the choice of reference has a negligible effect.\n", "regional_smb_weights = models[\"carra_14\"].sel(\n", " rgi_id=ds_smb.rgi_id, time=ds_smb.time_int).area.fillna(0)\n", "regional_smb_obs = ds_smb.weighted(regional_smb_weights).mean(dim=\"rgi_id\")[\"ba\"] * 1000\n", "\n", "\n", "def plot_regional_smb(ax, keys):\n", " for key in keys:\n", " regional_smb(models[key]).plot(ax=ax, label=CONFIGS[key])\n", " ax.plot(ds_smb.time_int, regional_smb_obs, c=\"k\", label=\"Regionally aggregated obs.\")\n", " ax.set_ylabel(\"Specific MB (mm w.e. yr$^{-1}$)\")\n", " ax.grid(True)\n", " ax.legend(bbox_to_anchor=[1.01, 0.5], loc=\"center left\")\n", "\n", "\n", "def smb_stats(model_series, obs_series):\n", " \"\"\"Bias, RMSD and correlation of a modelled SMB series against the SMB-map obs.\"\"\"\n", " aligned = model_series.sel(time=ds_smb.time_int)\n", " return {\n", " \"bias\": float(np.mean(aligned) - np.mean(obs_series)),\n", " \"rmsd\": float(np.sqrt(np.mean((aligned.values - obs_series.values) ** 2))),\n", " \"corr\": float(np.corrcoef(aligned.values, obs_series.values)[0, 1]),\n", " }\n", "\n", "\n", "def plot_single_glacier_smb(rgi_id, keys):\n", " fig, ax = plt.subplots()\n", " for key in keys:\n", " ds = models[key].sel(rgi_id=rgi_id)\n", " smb = ds.volume.diff(dim=\"time\", label=\"lower\").reindex(time=ds.time) / ds.area * 900\n", " smb.plot(ax=ax, label=CONFIGS[key])\n", " ax.plot(ds_smb.time_int, ds_smb[\"ba\"].sel(rgi_id=rgi_id) * 1000, c=\"k\", label=\"Observations\")\n", " ax.set_title(str(rgi_id.values))\n", " ax.set_ylabel(\"Specific MB (mm w.e. yr$^{-1}$)\")\n", " ax.grid(True)\n", " ax.legend()\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "f1c858f7", "metadata": {}, "source": [ "## Decision summary\n", "\n", "For each shortlisted configuration: the 2010 area mismatch against the outline\n", "observations (the independent check — 2025 is matched by construction, see note above,\n", "so it is not informative here), the 2000-2020 volume-change mismatch against Hugonnet et\n", "al. 2021 (adapted to the 2025 outlines, see `Landsvirkjun consultancy.md`), and the regional\n", "SMB fit against Andri's SMB maps. `hugonnet_mismatch_pct` is positive when the model\n", "loses *more* mass than observed, negative when it loses *less*." ] }, { "cell_type": "code", "execution_count": 4, "id": "452d673a", "metadata": { "execution": { "iopub.execute_input": "2026-07-07T21:44:06.779963Z", "iopub.status.busy": "2026-07-07T21:44:06.779824Z", "iopub.status.idle": "2026-07-07T21:44:06.885227Z", "shell.execute_reply": "2026-07-07T21:44:06.884380Z" } }, "outputs": [ { "data": { "text/html": [ "
| \n", " | area_2010_km2 | \n", "area_2010_mismatch_pct | \n", "mass_change_2000_2020_Gt | \n", "hugonnet_mismatch_pct | \n", "regional_smb_bias_mmwe | \n", "regional_smb_rmsd_mmwe | \n", "regional_smb_corr | \n", "
|---|---|---|---|---|---|---|---|
| config | \n", "\n", " | \n", " | \n", " | \n", " | \n", " | \n", " | \n", " |
| CARRA, trapez. λ=14 (~8°) | \n", "10427.07 | \n", "-1.73 | \n", "-175.24 | \n", "-6.39 | \n", "-375.16 | \n", "628.73 | \n", "0.41 | \n", "
| CARRA, trapez. λ=14 (~8°), Monthly TI | \n", "10470.89 | \n", "-1.32 | \n", "-193.45 | \n", "3.34 | \n", "-480.59 | \n", "671.30 | \n", "0.62 | \n", "
| ERA5, trapez. λ=8 (~14°) | \n", "10409.72 | \n", "-1.90 | \n", "-204.37 | \n", "9.17 | \n", "-319.75 | \n", "700.56 | \n", "0.24 | \n", "
| ERA5, trapez. λ=8 (~14°), Monthly TI | \n", "10408.16 | \n", "-1.91 | \n", "-196.03 | \n", "4.72 | \n", "-398.21 | \n", "698.14 | \n", "0.51 | \n", "