Creating DTC Glaciers EO-DT-Enhanced Data Cubes (L2)

Creating DTC Glaciers EO-DT-Enhanced Data Cubes (L2)#

For this example, we focus on the Brúarjökull outlet glacier of the Vatnajökull ice cap in Iceland.
Brúarjökull drains into catchments that feed Iceland’s largest hydropower system, making it an excellent case study for exploring how variations in glacier mass balance influence seasonal and long-term water availability (see this interactive visualisation story on the DestinE website).

Conceptually, Level 2 (L2) data cubes build upon the observational foundation provided by L1 data cubes by integrating EO observations with physically based glacier modelling and data assimilation.
Whilst L1 data cubes only contain harmonised observations, L2 data cubes provide model-consistent estimates of glacier state and fluxes, together with associated uncertainty information.

In this notebook, we use Brúarjökull to demonstrate how L1 observational data cubes (created in this notebook) are ingested into the DTCG API data assimilation pipeline to generate L2 data cubes.
This process illustrates how new EO observations can be systematically combined with modelling frameworks to produce temporally continuous and physically consistent representations of glacier evolution.

The example highlights how L2 data cubes bridge the gap between observations and projections, and provide a basis for improved assessments of future glacier change and its implications for water availability in glacier-influenced hydropower systems.

If required, install the DTCG API using the following command:

!pip install 'dtcg[jupyter] @ git+https://github.com/DTC-Glaciers/dtcg'

Run this command in a notebook cell.

# Imports
import xarray as xr
import matplotlib.pyplot as plt

import dtcg.integration.oggm_bindings as oggm_bindings
import dtcg.integration.calibration as calibration
from dtcg import DEFAULT_L1_DATACUBE_URL, DEFAULT_L2_DATACUBE_URL
from dtcg.api.external.call import StreamDatacube

Exploring L2 data cubes#

IMPORTANT: Please note that the L2 data cubes generated in this notebook use a smaller ensemble size for the Monte Carlo simulation. Therefore, the resulting uncertainties are not comparable to those of the preprocessed L2 data cubes. For this reason, we load the equivalent preprocessed L2 data cubes generated above. Uncomment the next cell only if you really want to inspect the L2 data cubes generated in this notebook.

def get_l2_data_tree(rgi_id):
    return xr.open_datatree(
        f"{DEFAULT_L2_DATACUBE_URL}{rgi_id}.zarr",
        chunks={},
        engine="zarr",
        consolidated=True,
        decode_cf=True,
    )


datacube_handler = GeoZarrHandler(ds=get_l2_data_tree(rgi_id_ice))

After generating the L2 data cubes and adding them to the data tree, we can now explore the results.
In general, each L2 data cube contains three datasets:

'monthly'
'annual_hydro'
'daily_smb'

They contain different variables at different temporal resolutions.

Let’s look at the monthly glacier volume evolution for our two realisations:

datacube_handler.data_tree["L2_Daily_Hugonnet_2010_2020"]["monthly"]

<xarray.DataTree 'monthly'>
Group: /L2_Daily_Hugonnet_2010_2020/monthly
    Dimensions:                   (member: 8, time: 310, rgi_id: 1)
    Coordinates:
      * member                    (member) <U7 224B 'Control' '0.05' ... '0.95'
      * time                      (time) float64 2kB 2e+03 2e+03 ... 2.026e+03
        calendar_month            (time) int64 2kB dask.array<chunksize=(310,), meta=np.ndarray>
        calendar_year             (time) int64 2kB dask.array<chunksize=(310,), meta=np.ndarray>
        hydro_month               (time) int64 2kB dask.array<chunksize=(310,), meta=np.ndarray>
        hydro_year                (time) int64 2kB dask.array<chunksize=(310,), meta=np.ndarray>
      * rgi_id                    (rgi_id) <U14 56B 'RGI60-06.00377'
    Data variables:
        area                      (member, time, rgi_id) float64 20kB dask.array<chunksize=(1, 310, 1), meta=np.ndarray>
        length                    (member, time, rgi_id) float64 20kB dask.array<chunksize=(1, 310, 1), meta=np.ndarray>
        specific_mb               (member, time, rgi_id) float64 20kB dask.array<chunksize=(1, 310, 1), meta=np.ndarray>
        specific_mb_calendar_cum  (member, time, rgi_id) float64 20kB dask.array<chunksize=(1, 310, 1), meta=np.ndarray>
        volume                    (member, time, rgi_id) float64 20kB dask.array<chunksize=(1, 310, 1), meta=np.ndarray>
    Attributes:
        Conventions:           CF-1.12
        comment:               The DTC Glaciers project is developed under the Eu...
        date_created:          2026-01-11T14:01:26.766555
        RGI-ID:                RGI60-06.00377
        glacier_attributes:    {}
        title:                 Datacube of observation-informed modelled variables.
        summary:               Observation-informed modelled variables for RGI6-I...
        calibration_strategy:  OGGM mass-balance model 'DailyTIModel' calibrated ...

import matplotlib.pyplot as plt
import seaborn as sns


def plot_quantile_bands(
    ds,
    x_dim,
    y_var,
    y_scaling,
    ax=None,
    color="C0",
    alpha=(0.3, 0.25, 0.2),
    label="median",
):

    x = ds[x_dim]
    # plot quantile bands
    bands = [
        (0.25, 0.75, alpha[0]),
        (0.15, 0.85, alpha[1]),
        (0.05, 0.95, alpha[2]),
    ]
    for q_low, q_high, a in bands:
        ax.fill_between(
            x,
            ds[y_var].sel(member=str(q_low)).isel(rgi_id=0) * y_scaling,
            ds[y_var].sel(member=str(q_high)).isel(rgi_id=0) * y_scaling,
            color=color,
            alpha=a,
            linewidth=0,
        )

    # Median line
    ax.plot(
        x,
        ds[y_var].sel(member="0.5") * y_scaling,
        color=color,
        label=label,
    )

    return ax

colors = sns.color_palette("colorblind")

fig, ax = plt.subplots()
for L2_name, c in zip(
    ["L2_Daily_Hugonnet_2010_2020", "L2_Daily_Cryosat_2011_2020"], colors
):
    plot_quantile_bands(
        datacube_handler.data_tree[L2_name]["monthly"],
        "time",
        "volume",
        1e-9,
        ax=ax,
        label=L2_name,
        color=c,
    )
ax.legend()
ax.grid("on")
ax.set_ylabel("Volume (km³)")
ax.set_xlabel("Year")
ax.set_title("Monthly L2 data cube example with uncertainties")
plt.show()

../_images/35898e11c2a2c6408ee2cf42ff2347fde3c7fa0858f7404b2f940f6062436179.png

The shaded areas show the quantile bands (0.25–0.75), (0.15–0.85), and (0.05–0.95), which represent 50%, 70%, and 90% of the ensemble outputs from the MCS.

We can also examine glacier runoff at annual or monthly resolution:

fig, ax = plt.subplots()
for L2_name, c in zip(
    ["L2_Daily_Hugonnet_2010_2020", "L2_Daily_Cryosat_2011_2020"], colors
):
    plot_quantile_bands(
        datacube_handler.data_tree[L2_name]["annual_hydro"],
        "time",
        "runoff",
        1e-9,
        ax=ax,
        label=L2_name,
        color=c,
    )

ax.legend()
ax.grid("on")
ax.set_ylabel("Runoff (Mt)")
ax.set_xlabel("Year")
ax.set_title("Annual L2 data cube example with uncertainties")
plt.show()

../_images/d4550a827cce28a6cffd075b5417f14f1adb2b41cf481e4d4b30c9c4b92d414a.png

fig, axs = plt.subplots(1, 2, figsize=(10, 5))

for ax, L2_name in zip(
    axs, ["L2_Daily_Hugonnet_2010_2020", "L2_Daily_Cryosat_2011_2020"]
):
    monthly_runoff = datacube_handler.data_tree[L2_name][
        "annual_hydro"
    ].runoff_monthly.sel(member="0.5")
    monthly_runoff = (
        monthly_runoff.rolling(time=31, center=True, min_periods=1).mean() * 1e-9
    )
    monthly_runoff.clip(0).plot(
        ax=ax,
        cmap="Blues",
        cbar_kwargs={
            "label": "Runoff (Mt)",
        },
        vmax=1000,
    )
    ax.set_title(L2_name)
    ax.set_xlabel("Months")
    ax.set_ylabel("Years")
plt.suptitle("Monthly hydro L2 data cube example, only median shown")
plt.tight_layout()

../_images/5323b30155940f0472358a2e47ac8f3125b56416c0d3c9858d99946cfc6098a7.png

Precomputed L2 data cubes#

We prepared preprocessed L2 data cubes for our case study regions in Iceland and Austria (see this notebook for an overview of the available glaciers).

These preprocessed L2 data cubes, together with the corresponding L1 data cubes, are stored as GeoZarr and can be accessed either by streaming directly with the DTCG API by providing an RGI6 ID, or the name of a glacier:

from dtcg import DEFAULT_L2_DATACUBE_URL
from dtcg.api.external.call import StreamDatacube

streamer = StreamDatacube(server=DEFAULT_L2_DATACUBE_URL)
datacube = streamer.stream_datacube(glacier=rgi_id_ice)

or using the following function:

from dtcg import DEFAULT_L2_DATACUBE_URL


def get_l2_data_tree(rgi_id):
    return xr.open_datatree(
        f"{DEFAULT_L2_DATACUBE_URL}{rgi_id}.zarr",
        chunks={},
        engine="zarr",
        consolidated=True,
        decode_cf=True,
    )

Creating DTC Glaciers EO-DT-Enhanced Data Cubes (L2)

Contents

Creating DTC Glaciers EO-DT-Enhanced Data Cubes (L2)#

Assimilating EO-data within DTC-Glaciers#

Preparing L1 data for assimilation#

Steps#

Climate forcing data from DestinE Platform#

Data assimilation and uncertainty propagation with Monte Carlo Simulation#

Running the assimilation and generating L2 data cubes#

Adding L2 data cubes to GeoZarr#

Exploring L2 data cubes#

Precomputed L2 data cubes#