Glaciers as water resources: part 1 (idealized climate)

Glaciers as water resources: part 1 (idealized climate)#

Our previous notebooks showed how to use OGGM to answer theoretical questions about glaciers, using idealized experiments: the effect of a different slope, the concept of the equilibrium line altitude, the mass balance gradient, etc. Now, how do we use OGGM to explore real glaciers? This notebook gives us some insight.

Goals of this notebook:

prepare a model run for a real world glacier
run simulations using different idealized climate scenarios to explore the role of glaciers as water resources
understand the concept of “peak water”

This version of the notebook works for OGGM versions 1.6 and above. If you are encountering errors, you are probably running an earlier version. The earlier version of this notebook can be found at: glacier_water_resources_pre16.ipynb

Setting the scene: glacier runoff and “peak water”#

“If glaciers melt, there won’t be water in mountain anymore”.

This is a sentence that we hear often from people we meet, or sometimes even in news articles. In fact, the role of glaciers in the hydrological cycle is more complex than that. In this notebook, we will explore this question using idealized climate scenarios applied to real glaciers!

Before we continue, let’s have a look at the expected contribution of glaciers to local streamflow. The chart below shows an idealized scenario where the climate is first constant (t\(_0\)-t\(_1\), and the glacier is in equilibrium with its climate) and then a warming occurs (t\(_1\)-t\(_2\)). This graph makes a few very important points, which we will explore together in this notebook.

When a glacier is in equilibrium, a glacier does not contribute to the annual runoff at all.
When the climate is warming, glaciers are losing mass. This water contributes to downstream runoff, and the runoff increases.
If climate warms even more, glaciers will continue to lose mass and become significantly smaller. When there isn’t much ice left to melt each year (or when climate stabilizes), their contribution will become smaller until becoming zero again.
In the new equilibrium, the annual runoff is the same as before, but the seasonal contribution changed.

We will now get back to all these points together, using OGGM!

Fig 1 from https://www.nature.com/articles/s41558-017-0049-x

Graphic from Huss & Hock (2018)

Setup#

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_context('notebook')  # plot defaults

import xarray as xr
import salem
import numpy as np
import pandas as pd

import oggm.cfg
from oggm import utils, workflow, tasks, graphics
from oggm.sandbox.edu import run_constant_climate_with_bias

# OGGM options
oggm.cfg.initialize(logging_level='WARNING')
oggm.cfg.PATHS['working_dir'] = utils.gettempdir(dirname='WaterResources')
oggm.cfg.PARAMS['min_ice_thick_for_length'] = 1  # a glacier is defined when ice is thicker than 1m
oggm.cfg.PARAMS['store_model_geometry'] = True

2023-08-27 22:01:41: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-08-27 22:01:41: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-08-27 22:01:41: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-08-27 22:01:42: oggm.cfg: PARAMS['min_ice_thick_for_length'] changed from `0.0` to `1`.
2023-08-27 22:01:42: oggm.cfg: PARAMS['store_model_geometry'] changed from `False` to `True`.

Define the glacier we will play with#

For this notebook we use the Hintereisferner, Austria. Some other possibilities to play with:

Hintereisferner, Austria: RGI60-11.00897
Artesonraju, Peru: RGI60-16.02444
Rikha Samba, Nepal: RGI60-15.04847
Parlung No. 94, China: RGI60-15.11693

And virtually any glacier you can find the RGI Id from, e.g. in the GLIMS viewer! Large glaciers may need longer simulations to see changes though. For less uncertain calibration parameters, we also recommend to pick one of the many reference glaciers in this list, where we make sure that observations of mass balance are better matched.

Let’s start with Hintereisferner first and you’ll be invited to try out your favorite glacier at the end of this notebook.

# Hintereisferner
rgi_id = 'RGI60-11.00897'

Preparing the glacier data#

This can take up to a few minutes on the first call because of the download of the required data:

# We pick the elevation-bands glaciers because they run a bit faster - but they create more step changes in the area outputs
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/elev_bands/W5E5'
gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=5, prepro_border=80, prepro_base_url=base_url)[0]

2023-08-27 22:01:43: oggm.workflow: init_glacier_directories from prepro level 5 on 1 glaciers.
2023-08-27 22:01:43: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

Interactive glacier map#

A first glimpse on the glacier of interest.

Tip: You can use the mouse to pan and zoom in the map

# One interactive plot below requires Bokeh
try:
    import holoviews as hv
    hv.extension('bokeh')
    import geoviews as gv
    import geoviews.tile_sources as gts

    sh = salem.transform_geopandas(gdir.read_shapefile('outlines'))
    out = (gv.Polygons(sh).opts(fill_color=None, color_index=None) * 
           gts.tile_sources['EsriImagery'] * 
           gts.tile_sources['StamenLabels']).opts(width=800, height=500, active_tools=['pan', 'wheel_zoom'])
except:
    # The rest of the notebook works without this dependency
    out = None
out

For OGGM, glaciers are “1.5” dimensional along their flowline:

fls = gdir.read_pickle('model_flowlines')
graphics.plot_modeloutput_section(fls);

../_images/1a8fc939201eb7b14bbfc93af2e06a4998b4c658676f295de560842ce1fdb2e9.png

Climate change simulation#

We start by generating a temperature change scenario. We define a short period of constant temperature (at -2°C, to be consistent with above), followed by a linear increase (2.5 degrees in about 150 years) and a stabilization:

ny_s = 50  # Start
ny_t = 150  # Trend
ny_e = 102  # Stabilisation
temp_bias_ts = np.concatenate([np.full(ny_s, -2.), np.linspace(-2, 0.5, ny_t), np.full(ny_e, 0.5)])
temp_bias_ts = pd.Series(temp_bias_ts, index=np.arange(ny_s + ny_t + ny_e))
temp_bias_ts.plot(); plt.xlabel('Year'); plt.ylabel('Temperature bias (°C)');

../_images/1157deec75da40fe8f28456a70dba85fe64ee069549bb7365beba1283919b9b1.png

# file identifier where the model output is saved
file_id = '_lin_temp'

# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs
tasks.run_with_hydro(gdir,
                     temp_bias_ts=temp_bias_ts,  # the temperature bias timeseries we just created
                     run_task=run_constant_climate_with_bias,  # which climate scenario? See following notebook for other examples
                     y0=2009, halfsize=10,  # Period which we will average and constantly repeat
                     store_monthly_hydro=True,  # Monthly ouptuts provide additional information
                     init_model_filesuffix='_spinup',  # We want to start from the glacier in equibrium we created earlier
                     output_filesuffix=file_id);  # an identifier for the output file, to read it later

Let’s have a look at our glacier evolution:

with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix=file_id)) as ds:
    # The last step of hydrological output is NaN (we can't compute it for this year)
    ds = ds.isel(time=slice(0, -1)).load()

fig, axs = plt.subplots(nrows=3, figsize=(10, 7), sharex=True)
ds.volume_m3.plot(ax=axs[0]);
ds.length_m.plot(ax=axs[1]);
temp_bias_ts.plot(ax=axs[2], c='C3');
axs[0].set_xlabel(''); axs[0].set_title(f'{rgi_id}'); axs[1].set_xlabel(''); axs[2].set_xlabel('Years'); axs[2].set_ylabel('Temperature bias (°C)');

../_images/d76f863c9e21ad0d7a72943a5c405880eedd9b738d03f4299e891cf0570c7e4b.png

In this scenario, the glacier would melt almost entirely (it looses about 90% of its volume). Note also that the glacier continues to adjust after the temperature has stabilized. What are the implications for downstream runoff?

Annual runoff#

As glaciers retreat, they contribute to sea level rise (visit the World Glaciers Explorer OGGM-Edu app for more information!). This is not what we are interested in here. Indeed, they will also have important local impacts: in this notebook, we will have a look at their impact on streamflow.

Let’s take a look at some of the hydrological outputs computed by OGGM. We start by creating a pandas DataFrame of all “1D” (annual) variables in the output dataset:

sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]
df_annual = ds[sel_vars].to_dataframe()

Then we can select the hydrological variables and sum them to get the total annual runoff:

# Select only the runoff variables
runoff_vars = ['melt_off_glacier', 'melt_on_glacier', 'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
# Convert them to megatonnes (instead of kg)
df_runoff = df_annual[runoff_vars] * 1e-9
# We smooth the output, which is otherwize noisy because of area discretization
df_runoff = df_runoff.rolling(31, center=True, min_periods=1).mean()
fig, ax = plt.subplots(figsize=(10, 3.5), sharex=True)
df_runoff.sum(axis=1).plot(ax=ax);
plt.ylabel('Runoff (Mt)'); plt.xlabel('Years'); plt.title(f'Total annual runoff for {rgi_id}');

../_images/faae4c7e9006ee7cb94aaee694c857d1eb75b3b4435ee951f89aaeadd65aa801.png

The hydrological variables are computed on the largest possible area that was covered by glacier ice during the simulation. This is equivalent to the runoff that would be measured at a fixed-gauge hydrological station at the glacier terminus.

The total annual runoff consists of the following components:

melt off-glacier: snow melt on areas that are now glacier free (i.e. 0 in the year of largest glacier extent, in this example at the start of the simulation)
melt on-glacier: ice + seasonal snow melt on the glacier
liquid precipitaton on- and off-glacier (the latter being zero at the year of largest glacial extent, in this example at start of the simulation)

f, ax = plt.subplots(figsize=(10, 6));
df_runoff.plot.area(ax=ax, color=sns.color_palette('rocket')); plt.xlabel('Years'); plt.ylabel('Runoff (Mt)'); plt.title(rgi_id);

../_images/aae0c8875c36dc7fc8a325e1ae0fa5806a85ffffc21d6cc983fe9a48d24676f8.png

The glacier length and volume decrease from year 50 onwards until about year 250 - this is the glacier retreat phase. Afterwards, length and volume stabilize at a constant value indicating that the glacier has reached equilibrium.

Now study the graph above and it’s four main components.

Questions to address:

When is “peak water” reached? Does this have anything to do with a change in the temperature trend itself? If not, what is the reason for “peak water” to occur?
Verify that the total fixed-gauge runoff (remember what this means?) is the same at the beginning and the end of the simulation. Why is that so? What changed in between?
The temperature stabilizes at year 200, yet all hydrological variables continue to change after that. What is happening there?
What is the contribution of liquid precipitation at the beginning and the end of the simulation? What changed?
Can you verify that the net glacier contribution to runoff (i.e. - \(\Delta M / \Delta t\)) is zero at the beginning and the end of the simulation and positive in between? Why do we expect this behavior?

Hint: compare the results to the idealised runoff graphic that we introduced at the beginning of this chapter

# Solution of the last question above
glacier_mass = ds.volume_m3.to_series() * oggm.cfg.PARAMS['ice_density'] * 1e-9  # In Megatonnes, Mt
glacier_mass = glacier_mass.rolling(31, center=True, min_periods=1).mean()
(- glacier_mass.diff()).plot()
plt.axhline(y=0, color='k', ls=':')
plt.ylabel('Annual glacier mass change (Mt yr$^{-1}$)')
plt.xlabel('Years'); plt.title('Glacier contribution to annual runoff');

../_images/d76d699bc8264d517a51d2c97be2741040cc3a2e93d1a52b4c06b78c2aeb66d6.png

Monthly runoff#

The “2D” variables contain the same hydrological data, but at monthly resolution (dimensions [time, month]). For example, monthly runoff can be computed as:

# Select only the runoff variables and convert them to megatonnes (instead of kg)
monthly_runoff = ds['melt_off_glacier_monthly'] + ds['melt_on_glacier_monthly'] + ds['liq_prcp_off_glacier_monthly'] + ds['liq_prcp_on_glacier_monthly']
monthly_runoff = monthly_runoff.rolling(time=31, center=True, min_periods=1).mean() * 1e-9
monthly_runoff.clip(0).plot(cmap='Blues', cbar_kwargs={'label': 'Runoff (Mt)'}); plt.xlabel('Months'); plt.ylabel('Years');

../_images/8138a25f7dafde690d8c7428b31b205f7e7ea946b145b0f9f1aa06234b97b4c4.png

As we can see, the runoff is approximately zero during the winter months, while relatively high during the summer months.

Now let’s compare the actual runoff to total precipitation over the basin:

# Compute total precipitation (Snow + Liquid)
tot_precip = ds['liq_prcp_off_glacier_monthly'] + ds['liq_prcp_on_glacier_monthly'] + ds['snowfall_off_glacier_monthly'] + ds['snowfall_on_glacier_monthly']
tot_precip *= 1e-9  # in Mt

# Plot these data at year 0
yr = 0
r = monthly_runoff.sel(time=yr)
p = tot_precip.sel(time=yr)

f, ax = plt.subplots(figsize=(10, 6));
r.plot(ax=ax, color='C3', label='Monthly runoff', linewidth=3);
p.plot(ax=ax, color='C0', label='Monthly precipitation', linewidth=3);
ax.fill_between(r.month_2d, r, p, where=(p >= r), facecolor='C0', interpolate=True, alpha=0.5)
ax.fill_between(r.month_2d, r, p, where=(r > p), facecolor='C3', interpolate=True, alpha=0.5)
plt.ylabel('Mt yr$^{-1}$'); plt.legend(loc='best');
plt.xlabel('Month'); plt.title(f'Total monthly runoff and precipitation at year {yr}');

../_images/8f4e4f67cae7661a269d33c0cf72843c65288f3587ec74e460440aa5b07a9ed1.png

At the begining of the simulation, the glacier is in equilibrium and its annual contribution to runoff is zero. In fact, this means that the blue area in the plot above is equal to the red area: the glaciated basin releases water during the summer, that accumulated in form of snow in the winter months.

In this region (Ötztal Alps, Austria), precipitation is relatively constant throughout the year. Discuss what the implications could be in other climates of the world, in particular in climates with a strong seasonality of precipitation.

Now let’s have a look at the seasonal change in runoff with time:

f, ax = plt.subplots(figsize=(10, 6));
cmap = sns.color_palette('magma', 3)
for i, yr in enumerate([0, 120, 300]):
    monthly_runoff.sel(time=yr).plot(ax=ax, color=cmap[i], label=f'Year {yr}')
plt.ylabel('Mt yr$^{-1}$'); plt.legend(loc='best');
plt.xlabel('Month'); plt.title('Total monthly runoff change with time');

../_images/7408aea3fb49891ed6a2684acf0db4e270caef0f333decaf19381ec69461fa2c.png

Questions to explore:

verify that the annual runoff is the same at year 0 than at year 300. What changed in between?
what is responsible for the switch from summer melt to spring melt?
how does the runoff change at “peak water” (approx. year 120)?
discuss the implications of this shift in seasonality for Hintereisferner and for other regions of the world.

Now repeat this notebook (ideally by copying it first) with other regions of the world, where climate might be very different.

Note that our explanations are very glacier-specific. So if you change the glacier by using another rgi_id, you also might need to adapt the explanations!

You can also change the parameters of the simulations!

What happens when the temperature change is larger/faster, smaller/slower? What are the implications for “peak water”?

Wrapping up#

Now let’s conclude this notebook by reproducing the famous plot by Huss & Hock (2018) with our data.

Graphic from Huss & Hock (2018)

# Make a dataframe out of the xarray dataset for simplicity
df_monthly_runoff = pd.DataFrame(monthly_runoff.data, index=ds.time, columns=monthly_runoff.month_2d)

# Create new columns for each season (we define seasons from the perspective of the northern Hemisphere)
df_monthly_runoff['Spring (Apr-Jun)'] = df_monthly_runoff[[4, 5, 6]].sum(axis=1)
df_monthly_runoff['Summer (Jul-Sep)'] = df_monthly_runoff[[7, 8, 9]].sum(axis=1)
df_monthly_runoff['Winter (Oct-Mar)'] = df_monthly_runoff[[10, 11, 12, 1, 2, 3]].sum(axis=1)
df_monthly_runoff.columns.name = 'Season'

# Plot it
# Attention the seasons and labels might need to be adapted if a glacier in a different region is used!
f, ax = plt.subplots(figsize=(10, 6));
df_monthly_runoff[['Spring (Apr-Jun)', 'Summer (Jul-Sep)', 'Winter (Oct-Mar)']].plot.area(ax=ax, color=sns.color_palette('rocket'));
plt.xlabel('Years'); plt.ylabel('Runoff (Mt)'); plt.title('Runoff by season');

../_images/d04f825a04fc6708a756e552fb9151775c158bfbb7d6fab72a0be9cbbf25f400.png

And now the actual plot:

f, (ax1, ax2) = plt.subplots(2, 1, gridspec_kw={'height_ratios': [2, 1]}, figsize=(10, 8), sharex=True)

p1 = df_monthly_runoff[['Spring (Apr-Jun)', 'Summer (Jul-Sep)', 'Winter (Oct-Mar)']].sum(axis=1)
p2 = df_monthly_runoff['Summer (Jul-Sep)'] # July to September is here the melt season
p1.plot(ax=ax1, color='C0', label='Annual runoff', linewidth=3)
p2.plot(ax=ax1, color='C3', label='Melt-season runoff', linewidth=3);
ax1.hlines([p1.loc[0], p2.loc[0]], 0, 300, color=['C0', 'C3'], linestyle=[':']);
ax1.legend(); ax1.set_xlabel(''); ax1.set_ylabel('Runoff (Mt)');

(temp_bias_ts + 2).plot(ax=ax2, color='C2', label='Temperature', linewidth=3);
ax2.legend(loc='lower right'); ax2.set_xlabel('Years'); ax2.set_ylabel('$\Delta T (°C)$');
sns.despine();

../_images/9e9c2f21e072694b7949ded8bb161a0402c9ff07000b787d9bd9f120379cba36.png

Note that the two plots above are very glacier-specific. If you use a glacier with another climate, you might need to adapt the plots, because we are always interested in those months where the runoff from the glacier contributes most to the total runoff. Example given, for Artesonraju in Peru, we are most interested in June, July, August and September as this is the dry season, where it almost does not rain at all.

Take home points#

Glaciers in equilibrium are not net water resources: they gain as much mass as they release
However, they have a seasonal buffer role, releasing water during the melt months
When glaciers melt, they become net water resources. “Peak water” is the point in time when glacier melt supply reaches its maximum, i.e. when the maximum runoff occurs.

References#

Huss, M. and Hock, R.: Global-scale hydrological response to future glacier mass loss, Nat. Clim. Chang., 8(2), 135–140, doi:10.1038/s41558-017-0049-x, 2018.
A recent blog post about the role of glaciers in the hydrological cycle

What’s next?#

Back to the table of contents