Hydrological mass-balance output#

New in version 1.5!

In two recent PRs (GH1224, GH1232 and GH1242), we have added a new task in OGGM, run_with_hydro, which adds mass-balance and runoff diagnostics to the OGGM output files.

This task is still experimental - it is tested for consistency and mass-conservation and we trust its output, but its API and functionalites might change in the future (in particular to make it faster and to add more functionality, as explained below).

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

# Make pretty plots
sns.set_style('ticks')
sns.set_context('notebook')

from oggm import cfg, utils, workflow, tasks, graphics

cfg.initialize(logging_level='WARNING')
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGMHydro')
cfg.PARAMS['store_model_geometry'] = True  # Hydro outputs needs the full model geometry

2023-03-07 12:27:28: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-03-07 12:27:28: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-07 12:27:28: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-03-07 12:27:28: 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:

Shallap Glacier: RGI60-16.02207
Artesonraju: RGI60-16.02444 (reference glacier)
Hintereisferner: RGI60-11.00897 (reference glacier)

And virtually any glacier you can find the RGI Id from, e.g. in the GLIMS viewer.

# 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
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc3/pcp2.5/no_match'
gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=5, prepro_border=80, prepro_base_url=base_url)[0]

2023-03-07 12:27:28: oggm.workflow: init_glacier_directories from prepro level 5 on 1 glaciers.
2023-03-07 12:27:28: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

Random climate commitment run#

Same as before, but this time randomly shuffling each year of the last 11 years:

# file identifier where the model output is saved
file_id = '_rd'
tasks.run_with_hydro(gdir, run_task=tasks.run_random_climate, nyears=100, y0=2014, halfsize=5, store_monthly_hydro=True,
                     seed=0, unique_samples=True,
                     output_filesuffix=file_id);

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

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()

Annual runoff#

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

# Select only the runoff variables and convert them to megatonnes (instead of kg)
runoff_vars = ['melt_off_glacier', 'melt_on_glacier', 'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
df_runoff = df_annual[runoff_vars] * 1e-9

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);

Monthly runoff#

ds_roll = ds.roll(month_2d=ds['calendar_month_2d'].data[0]-1, roll_coords=True)
ds_roll['month_2d'] = ds_roll['calendar_month_2d']

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

# Pick the variables we need (the 2d ones)
sel_vars = [v for v in ds_roll.variables if 'month_2d' in ds_roll[v].dims]

# Pick the first decade and average it
df_m_s = ds_roll[sel_vars].isel(time=slice(0, 10)).mean(dim='time').to_dataframe() * 1e-9
# Rename the columns for readability
df_m_s.columns = [c.replace('_monthly', '') for c in df_m_s.columns]
# Because of floating point precision sometimes runoff can be slightly below zero, clip
df_m_s = df_m_s.clip(0)

# Same for end
df_m_e = ds_roll[sel_vars].isel(time=slice(-11, -1)).mean(dim='time').to_dataframe() * 1e-9
df_m_e.columns = [c.replace('_monthly', '') for c in df_m_s.columns]
df_m_e = df_m_e.clip(0)

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 7), sharey=True);
df_m_s[runoff_vars].plot.area(ax=ax1, legend=False, title='Year 0-10', color=sns.color_palette("rocket"));
df_m_e[runoff_vars].plot.area(ax=ax2, title='Year 90-100', color=sns.color_palette("rocket"));
ax1.set_ylabel('Monthly runoff (Mt)'); ax1.set_xlabel('Month'); ax2.set_xlabel('Month');

CMIP5 projection run#

This time, let’s start from the estimated glacier state in 2020 and run a projection from there. See the run_with_gcm tutorial for more info.

For illustrative purposes, we also show how you can use two new keywords added in OGGM v1.5.3 to:

compute hydrological outputs with a fixed geometry prior to the RGI date
compute the hydrological outputs with the same reference glacier area in the projections and the past simulation

For an (experimental) use of a dynamic spinup instead of this fixed geometry spinup, visit dynamical_spinup.ipynb.

“Historical” run for the reference climate period#

# Run the historical period with a fixed geometry spinup
file_id = '_hist_hydro'
tasks.run_with_hydro(gdir, run_task=tasks.run_from_climate_data,
                     fixed_geometry_spinup_yr=1990,  # this tells to do spinup
                     ref_area_from_y0=True,  # 
                     store_monthly_hydro=True,
                     output_filesuffix=file_id);

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

# Read the data
with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix=file_id)) as ds:
    ds = ds.isel(time=slice(0, -1)).load()
    sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]
    df_annual = ds[sel_vars].to_dataframe()

f, ax = plt.subplots(figsize=(9, 6))
df_annual[['volume_m3', 'area_m2']].plot(ax=ax, secondary_y='area_m2');
plt.suptitle(f'Glacier volume and area with fixed geometry spinup (RGI date: {gdir.rgi_date})');

The plot shows what OGGM has done: glacier area is left constant until the RGI date, after which the dynamical model is used. Glacier volume however does vary, and is consistent with the glacier surface mass-balance (albeit without geometry feedbacks). The hydrological output is also consistent with climate and a fixed glacier geometry:

# Select only the runoff variables and convert them to megatonnes (instead of kg)
f, ax = plt.subplots(figsize=(10, 6));
(df_annual[runoff_vars] * 1e-9).plot.area(ax=ax, color=sns.color_palette("rocket")); 
plt.xlabel('Years'); plt.ylabel('Runoff (Mt)'); plt.title(rgi_id);

Projection run#

# "Downscale" the climate data
from oggm.shop import gcm_climate
bp = 'https://cluster.klima.uni-bremen.de/~oggm/cmip5-ng/pr/pr_mon_CCSM4_{}_r1i1p1_g025.nc'
bt = 'https://cluster.klima.uni-bremen.de/~oggm/cmip5-ng/tas/tas_mon_CCSM4_{}_r1i1p1_g025.nc'
for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    # Download the files
    ft = utils.file_downloader(bt.format(rcp))
    fp = utils.file_downloader(bp.format(rcp))
    # bias correct them
    workflow.execute_entity_task(gcm_climate.process_cmip_data, [gdir],
                                 filesuffix='_CCSM4_{}'.format(rcp),  # recognize the climate file for later
                                 fpath_temp=ft,  # temperature projections
                                 fpath_precip=fp,  # precip projections
                                 );

2023-03-07 12:27:48: oggm.workflow: Execute entity tasks [process_cmip_data] on 1 glaciers
2023-03-07 12:28:02: oggm.workflow: Execute entity tasks [process_cmip_data] on 1 glaciers
2023-03-07 12:28:16: oggm.workflow: Execute entity tasks [process_cmip_data] on 1 glaciers
2023-03-07 12:28:30: oggm.workflow: Execute entity tasks [process_cmip_data] on 1 glaciers

for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    rid = '_CCSM4_{}'.format(rcp)
    tasks.run_with_hydro(gdir, run_task=tasks.run_from_climate_data,
                         climate_filename='gcm_data',  # use gcm_data, not climate_historical
                         climate_input_filesuffix=rid,  # use the chosen scenario
                         init_model_filesuffix='_hist_hydro',  # this is important! Start from 2020 glacier
                         ref_geometry_filesuffix='_hist_hydro',  # also use this as area reference
                         ref_area_from_y0=True,  # and keep the same reference area as for the historical simulations
                         output_filesuffix=rid,  # recognize the run for later
                         store_monthly_hydro=True,  # add monthly diagnostics
                         );

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

RCP2.6#

file_id = '_CCSM4_rcp26'

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()

ds.volume_m3.plot();

Annual runoff#

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

# Select only the runoff variables and convert them to megatonnes (instead of kg)
runoff_vars = ['melt_off_glacier', 'melt_on_glacier', 'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
df_runoff = df_annual[runoff_vars].clip(0) * 1e-9

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);

Note that “melt_off_glacier” is not zero at the beginning of the simulation - indeed the reference area is that of the 2003 glacier.

Monthly runoff#

ds_roll = ds.roll(month_2d=ds['calendar_month_2d'].data[0]-1, roll_coords=True)
ds_roll['month_2d'] = ds_roll['calendar_month_2d']

# Pick the variables we need (the 2d ones)
sel_vars = [v for v in ds_roll.variables if 'month_2d' in ds_roll[v].dims]

# Pick the first decade and average it
df_m_s = ds_roll[sel_vars].isel(time=slice(0, 10)).mean(dim='time').to_dataframe() * 1e-9
# Rename the columns for readability
df_m_s.columns = [c.replace('_monthly', '') for c in df_m_s.columns]
# Because of floating point precision sometimes runoff can be slightly below zero, clip
df_m_s = df_m_s.clip(0)

# Same for end
df_m_e = ds_roll[sel_vars].isel(time=slice(-11, -1)).mean(dim='time').to_dataframe() * 1e-9
df_m_e.columns = [c.replace('_monthly', '') for c in df_m_s.columns]
df_m_e = df_m_e.clip(0)

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 7), sharey=True);
df_m_s[runoff_vars].plot.area(ax=ax1, legend=False, title='Year 0-10', color=sns.color_palette("rocket"));
df_m_e[runoff_vars].plot.area(ax=ax2, title='Year 90-100', color=sns.color_palette("rocket"));
ax1.set_ylabel('Monthly runoff (Mt)'); ax1.set_xlabel('Month'); ax2.set_xlabel('Month');

RCP8.5#

file_id = '_CCSM4_rcp85'

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()

ds.volume_m3.plot();

Annual runoff#

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

# Select only the runoff variables and convert them to megatonnes (instead of kg)
runoff_vars = ['melt_off_glacier', 'melt_on_glacier', 'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
df_runoff = df_annual[runoff_vars] * 1e-9

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);

Monthly runoff#

ds_roll = ds.roll(month_2d=ds['calendar_month_2d'].data[0]-1, roll_coords=True)
ds_roll['month_2d'] = ds_roll['calendar_month_2d']

# Pick the variables we need (the 2d ones)
sel_vars = [v for v in ds_roll.variables if 'month_2d' in ds_roll[v].dims]

# Pick the first decade and average it
df_m_s = ds_roll[sel_vars].isel(time=slice(0, 10)).mean(dim='time').to_dataframe() * 1e-9
# Rename the columns for readability
df_m_s.columns = [c.replace('_monthly', '') for c in df_m_s.columns]
# Because of floating point precision sometimes runoff can be slightly below zero, clip
df_m_s = df_m_s.clip(0)

# Same for end
df_m_e = ds_roll[sel_vars].isel(time=slice(-11, -1)).mean(dim='time').to_dataframe() * 1e-9
df_m_e.columns = [c.replace('_monthly', '') for c in df_m_s.columns]
df_m_e = df_m_e.clip(0)

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 7), sharey=True);
df_m_s[runoff_vars].plot.area(ax=ax1, legend=False, title='Year 0-10', color=sns.color_palette("rocket"));
df_m_e[runoff_vars].plot.area(ax=ax2, title='Year 90-100', color=sns.color_palette("rocket"));
ax1.set_ylabel('Monthly runoff (Mt)'); ax1.set_xlabel('Month'); ax2.set_xlabel('Month');

Calculating peak water under different climate scenarios#

A typical usecase for simulating and analysing the hydrological ouptuts of glaciers are for peak water estimations. For instance, we might want to know the point in time when the annual total runoff from a glacier reaches its maximum under a certain climate scenario.

The total runoff is a sum of the melt and liquid precipitation, from both on and off the glacier. Peak water can be calculated from the 11-year moving average of the total runoff (Huss and Hock, 2018).

For Hintereisferner we have alredy run the simulations for the different climate scenarios, so we can sum up the runoff variables and plot the output.

# Create the figure
f, ax = plt.subplots(figsize=(18, 7), sharex=True)
# Loop all scenarios
for i, rcp in enumerate(['rcp26', 'rcp45', 'rcp60', 'rcp85']):
    file_id = f'_CCSM4_{rcp}'
    # Open the corresponding data.
    with xr.open_dataset(gdir.get_filepath('model_diagnostics',
                                           filesuffix=file_id)) as ds:
        # Load the data into a dataframe
        ds = ds.isel(time=slice(0, -1)).load()

    # Select annual variables
    sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]
    # And create a dataframe
    df_annual = ds[sel_vars].to_dataframe()

    # Select the variables relevant for runoff.
    runoff_vars = ['melt_off_glacier', 'melt_on_glacier',
                   'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
    # Convert to mega tonnes instead of kg.
    df_runoff = df_annual[runoff_vars].clip(0) * 1e-9
    # Sum the variables each year "axis=1", take the 11 year rolling mean
    # and plot it.
    running = df_runoff.sum(axis=1).rolling(window=11).mean()
    running.plot(ax=ax, label=rcp, color=sns.color_palette("rocket")[i])
ax.set_ylabel('Annual runoff (Mt)')
ax.set_xlabel('Year')
plt.title(rgi_id)
plt.legend();

For Hintereisferner, runoff continues to decrease throughout the 21st-century for all scenarios, indicating that peak water has already been reached sometime in the past. This is the case for many European glaciers. Let us take a look at another glacier, in a different climatical setting where a different hydrological projection can be expected.

We pick a glacier (RGI60-15.02420, unnamed) in the Eastern Himalayas. First we need to initialize its glacier directory.

# Unnamed glacier
rgi_id = 'RGI60-15.02420'

gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=5, prepro_border=80, prepro_base_url=base_url)[0]

2023-03-07 12:28:37: oggm.workflow: init_glacier_directories from prepro level 5 on 1 glaciers.
2023-03-07 12:28:37: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[51], line 1
----> 1 gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=5, prepro_border=80, prepro_base_url=base_url)[0]

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/workflow.py:533, in init_glacier_directories(rgidf, reset, force, from_prepro_level, prepro_border, prepro_rgi_version, prepro_base_url, from_tar, delete_tar)
    531     if cfg.PARAMS['dl_verify']:
    532         utils.get_dl_verify_data('cluster.klima.uni-bremen.de')
--> 533     gdirs = execute_entity_task(gdir_from_prepro, entities,
    534                                 from_prepro_level=from_prepro_level,
    535                                 prepro_border=prepro_border,
    536                                 prepro_rgi_version=prepro_rgi_version,
    537                                 base_url=prepro_base_url)
    538 else:
    539     # We can set the intersects file automatically here
    540     if (cfg.PARAMS['use_intersects'] and
    541             len(cfg.PARAMS['intersects_gdf']) == 0 and
    542             not from_tar):

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/workflow.py:191, in execute_entity_task(task, gdirs, **kwargs)
    187     if ng > 3:
    188         log.workflow('WARNING: you are trying to run an entity task on '
    189                      '%d glaciers with multiprocessing turned off. OGGM '
    190                      'will run faster with multiprocessing turned on.', ng)
--> 191     out = [pc(gdir) for gdir in gdirs]
    193 return out

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/workflow.py:191, in <listcomp>(.0)
    187     if ng > 3:
    188         log.workflow('WARNING: you are trying to run an entity task on '
    189                      '%d glaciers with multiprocessing turned off. OGGM '
    190                      'will run faster with multiprocessing turned on.', ng)
--> 191     out = [pc(gdir) for gdir in gdirs]
    193 return out

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/workflow.py:108, in _pickle_copier.__call__(self, arg)
    106 for func in self.call_func:
    107     func, kwargs = func
--> 108     res = self._call_internal(func, arg, kwargs)
    109 return res

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/workflow.py:102, in _pickle_copier._call_internal(self, call_func, gdir, kwargs)
     99     gdir, gdir_kwargs = gdir
    100     kwargs.update(gdir_kwargs)
--> 102 return call_func(gdir, **kwargs)

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/workflow.py:251, in gdir_from_prepro(entity, from_prepro_level, prepro_border, prepro_rgi_version, base_url)
    248 tar_base = utils.get_prepro_gdir(prepro_rgi_version, rid, prepro_border,
    249                                  from_prepro_level, base_url=base_url)
    250 from_tar = os.path.join(tar_base.replace('.tar', ''), rid + '.tar.gz')
--> 251 return oggm.GlacierDirectory(entity, from_tar=from_tar)

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2529, in GlacierDirectory.__init__(self, rgi_entity, base_dir, reset, from_tar, delete_tar)
   2525     raise RuntimeError('RGI Version not supported: '
   2526                        '{}'.format(self.rgi_version))
   2528 # remove spurious characters and trailing blanks
-> 2529 self.name = filter_rgi_name(name)
   2531 # region
   2532 reg_names, subreg_names = parse_rgi_meta(version=self.rgi_version[0])

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_funcs.py:734, in filter_rgi_name(name)
    728 def filter_rgi_name(name):
    729     """Remove spurious characters and trailing blanks from RGI glacier name.
    730 
    731     This seems to be unnecessary with RGI V6
    732     """
--> 734     if name is None or len(name) == 0:
    735         return ''
    737     if name[-1] in ['À', 'È', 'è', '\x9c', '3', 'Ð', '°', '¾',
    738                     '\r', '\x93', '¤', '0', '`', '/', 'C', '@',
    739                     'Å', '\x06', '\x10', '^', 'å', ';']:

TypeError: object of type 'numpy.float64' has no len()

Then, we have to process the climate data for the new glacier.

# Do we need to download data again? I guess it is stored somewhere,
# but this is an easy way to loop over it for bias correction.
bp = 'https://cluster.klima.uni-bremen.de/~oggm/cmip5-ng/pr/pr_mon_CCSM4_{}_r1i1p1_g025.nc'
bt = 'https://cluster.klima.uni-bremen.de/~oggm/cmip5-ng/tas/tas_mon_CCSM4_{}_r1i1p1_g025.nc'
for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    # Download the files
    ft = utils.file_downloader(bt.format(rcp))
    fp = utils.file_downloader(bp.format(rcp))
    workflow.execute_entity_task(gcm_climate.process_cmip_data, [gdir],
                                 # Name file to recognize it later
                                 filesuffix='_CCSM4_{}'.format(rcp),  
                                 # temperature projections
                                 fpath_temp=ft,
                                 # precip projections
                                 fpath_precip=fp,
                                 );

With this done, we can run the simulations for the different climate scenarios (for simplicity here, we don’t care about historical simulations).

for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    rid = '_CCSM4_{}'.format(rcp)
    tasks.run_with_hydro(gdir, run_task=tasks.run_from_climate_data, ys=2020,
                         # use gcm_data, not climate_historical
                         climate_filename='gcm_data',
                         # use the chosen scenario
                         climate_input_filesuffix=rid,
                         # this is important! Start from 2020 glacier
                         init_model_filesuffix='_historical',
                         # recognize the run for later
                         output_filesuffix=rid,
                         # add monthly diagnostics
                         store_monthly_hydro=True,
                        );

Now we can create the same plot as before in order to visualize peak water

# Create the figure
f, ax = plt.subplots(figsize=(18, 7))
# Loop all scenarios
for i, rcp in enumerate(['rcp26', 'rcp45', 'rcp60', 'rcp85']):
    file_id = f'_CCSM4_{rcp}'
    # Open the corresponding data in a context manager.
    with xr.open_dataset(gdir.get_filepath('model_diagnostics',
                                           filesuffix=file_id)) as ds:
        # Load the data into a dataframe
        ds = ds.isel(time=slice(0, -1)).load()

    # Select annual variables
    sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]
    # And create a dataframe
    df_annual = ds[sel_vars].to_dataframe()

    # Select the variables relevant for runoff.
    runoff_vars = ['melt_off_glacier', 'melt_on_glacier',
                   'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
    # Convert to mega tonnes instead of kg.
    df_runoff = df_annual[runoff_vars].clip(0) * 1e-9
    # Sum the variables each year "axis=1", take the 11 year rolling mean
    # and plot it.
    df_runoff.sum(axis=1).rolling(window=11).mean().plot(ax=ax, label=rcp,
                                  color=sns.color_palette("rocket")[i]
                                  )
ax.set_ylabel('Annual runoff (Mt)')
ax.set_xlabel('Year')
plt.title(rgi_id)
plt.legend();

Unlike for Hintereisferner, these projections indicate that the annual runoff will increase in all the scenarios for the first half of the century. The higher RCP scenarios can reach peak water later in the century, since the excess melt can continue to increase. For the lower RCP scenarios on the other hand, the glacier might be approaching a new equilibirum, which reduces the runoff earlier in the century (Rounce et. al., 2020). After peak water is reached (RCP2.6: ~2055, RCP8.5: ~2070 in these projections), the annual runoff begin to decrease. This decrease occurs because the shrinking glacier is no longer able to support the high levels of melt.

What’s next?#

return to the OGGM documentation
back to the table of contents

OGGM tutorials

Hydrological mass-balance output

Contents

Hydrological mass-balance output#

Define the glacier we will play with#

Preparing the glacier data#

“Commitment run”#

Annual runoff#

Monthly runoff#

Random climate commitment run#

Annual runoff#

Monthly runoff#

CMIP5 projection run#

“Historical” run for the reference climate period#

Projection run#

RCP2.6#

Annual runoff#

Monthly runoff#

RCP8.5#

Annual runoff#

Monthly runoff#

Calculating peak water under different climate scenarios#

What’s next?#