A quick look into the mass-balance calibration procedure

The default mass-balance model of OGGM is a very standard temperature index melt model. What’s … unusual about it, is the calibration procedure. We have explained it in quite some length in the documentation, in our paper, in the original publication that introduced it for the first time, and we have created a website to monitor its performance. The method is not perfect (by far), but it is quite powerful: I’ve often said that it is the best idea that Ben ever had, and he has a lot of good ideas.

However, experience shows that most people (including us sometimes ;) don’t understand how it works. Let’s add a new tutorial to the list, and jump on the opportunity to play around with the calibration procedure (don’t do this at home!).


import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import os

import oggm
from oggm import cfg, utils, workflow, tasks, graphics
from oggm.core import massbalance, climate
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-ref-mb', reset=True)
cfg.PARAMS['border'] = 10
2022-10-07 12:56:59: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2022-10-07 12:56:59: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2022-10-07 12:56:59: oggm.cfg: Multiprocessing: using all available processors (N=2)
2022-10-07 12:56:59: oggm.cfg: PARAMS['border'] changed from `40` to `10`.

We start from two well known glaciers in the Austrian Alps, Kesselwandferner and Hintereisferner:

gdirs = workflow.init_glacier_directories(['RGI60-11.00787', 'RGI60-11.00897'], from_prepro_level=3)
InvalidParamsError                        Traceback (most recent call last)
Cell In [3], line 1
----> 1 gdirs = workflow.init_glacier_directories(['RGI60-11.00787', 'RGI60-11.00897'], from_prepro_level=3)

File /usr/local/pyenv/versions/3.10.7/lib/python3.10/site-packages/oggm/workflow.py:351, in init_glacier_directories(rgidf, reset, force, from_prepro_level, prepro_border, prepro_rgi_version, prepro_base_url, from_tar, delete_tar)
    348     reset = utils.query_yes_no('Delete all glacier directories?')
    350 if from_prepro_level:
--> 351     url = utils.get_prepro_base_url(base_url=prepro_base_url,
    352                                     border=prepro_border,
    353                                     prepro_level=from_prepro_level,
    354                                     rgi_version=prepro_rgi_version)
    355     if cfg.PARAMS['has_internet'] and not utils.url_exists(url):
    356         raise InvalidParamsError("base url seems unreachable with these "
    357                                  "parameters: {}".format(url))

File /usr/local/pyenv/versions/3.10.7/lib/python3.10/site-packages/oggm/utils/_downloads.py:1278, in get_prepro_base_url(base_url, rgi_version, border, prepro_level)
   1275 """Extended base url where to find the desired gdirs."""
   1277 if base_url is None:
-> 1278     raise InvalidParamsError('Starting with v1.6, users now have to '
   1279                              'explicitly indicate the url they want '
   1280                              'to start from.')
   1282 if not base_url.endswith('/'):
   1283     base_url += '/'

InvalidParamsError: Starting with v1.6, users now have to explicitly indicate the url they want to start from.

The two glaciers are neighbors but have very different geometries:

f, ax = plt.subplots(figsize=(8, 8))
graphics.plot_googlemap(gdirs, ax=ax)

Calibration on glaciers with data

Like most models, they are always “perfect” where observations are available (both glaciers are reference glaciers in the WGMS database):

# Get HEF
gdir = gdirs[1]
# Get the reference mass-balance from the WGMS
ref_df = gdir.get_ref_mb_data()
# Get the calibrated mass-balance model
mbmod = massbalance.PastMassBalance(gdir)
# Compute the specific MB for this glacier
fls = gdir.read_pickle('inversion_flowlines')
ref_df['OGGM (calib)'] = mbmod.get_specific_mb(fls=fls, year=ref_df.index.values)
ref_df[['ANNUAL_BALANCE', 'OGGM (calib)']].plot();

Perfect? Not really: the model is calibrated so that the bias over the calibration period is zero, but other uncertainties remain: the variability in our timeseries is too high (this has something to do with precipitation amounts) and we have a stronger trend in the model than observations (this is due to the changing glacier geometry, that our model doesn’t know about here: so this is more a “feature” than a bug).

What’s important for OGGM here is the bias, which is zero:

ref_df[['ANNUAL_BALANCE', 'OGGM (calib)']].mean()

The calibration procedure did following: it found a year \(t^*\) in the past which is selected so that the glacier, with its present day geometry, would be in equilibrium with a 31-yr climate centered on \(t^*\) (see the documentation for more info). The year \(t^*\) and the associated temperature sensitivity \(\mu^*\) (units: mm K\(^{-1}\) yr\(^{-1}\)) for this glacier are:

gdir.read_json('local_mustar')['t_star'], gdir.read_json('local_mustar')['mu_star_glacierwide']

What about Kesselwandferner?

gdir = gdirs[0]
gdir.read_json('local_mustar')['t_star'], gdir.read_json('local_mustar')['mu_star_glacierwide']

These are very different values! Not surprising: the two glaciers are in the same climate (from the forcing data) but are very different in size, orientation and geometry. The resulting mass-balance values are also different:

ref_df = gdir.get_ref_mb_data()
mbmod = massbalance.PastMassBalance(gdir)
ref_df['OGGM (calib)'] = mbmod.get_specific_mb(fls=gdir.read_pickle('inversion_flowlines'), year=ref_df.index.values)
ref_df[['ANNUAL_BALANCE', 'OGGM (calib)']].mean()

Kesselwandferner has a mass-balance less negative than its neighbor, probably because it is smaller in size and spans less altitude.

Calibration on glaciers without data

Let’s trick OGGM a little. Let’s make it appear that no data is available at Kesselwandferner for calibration:

# Set an empty reference mass-balance

And re-run the mass-balance calibration:

cfg.PARAMS['run_mb_calibration'] = True
workflow.execute_entity_task(tasks.local_t_star, gdirs);
workflow.execute_entity_task(tasks.mu_star_calibration, gdirs);

Now, what’s the new \(\mu^*\) of Kesselwandferner?

gdir.read_json('local_mustar')['t_star'], gdir.read_json('local_mustar')['mu_star_glacierwide']

In the absence of other glaciers to interpolate from (usually the 10 closest), OGGM simply assigned the value of \(t^*\) from Hintereisferner to Kesselwandferner (1963). From this, \(\mu^*\) could be estimated. This “blind” mass-balance model is of course not as good as if we calibrated it using observations:

mbmod = massbalance.PastMassBalance(gdir)
ref_df['OGGM (blind)'] = mbmod.get_specific_mb(fls=gdir.read_pickle('inversion_flowlines'), year=ref_df.index.values)
ref_df[['ANNUAL_BALANCE', 'OGGM (calib)', 'OGGM (blind)']].mean()

Quite a negative bias! This is as good as OGGM can be in the absence of observations and without any other reference glacier in vicinity (OGGM needs at least one).

What else could we do to calibrate \(\mu^*\) for this glacier? We tried for example to simply interpolate \(\mu^*\), and found that on average, this leads to worse results. Here for example:

mbmod = massbalance.PastMassBalance(gdir, mu_star=221.275)  # Apply mu* from Hintereisferner
ref_df['OGGM (μ∗ from HEF)'] = mbmod.get_specific_mb(fls=gdir.read_pickle('inversion_flowlines'), year=ref_df.index.values)
ref_df[['ANNUAL_BALANCE', 'OGGM (calib)', 'OGGM (blind)', 'OGGM (μ∗ from HEF)']].mean()

Interpolating \(\mu^*\) instead of \(t^*\) is quite quite worse! And we found out that this is generally the case globally:

img Benefit of spatially interpolating \(t^*\) instead of \(\mu^*\) as shown by leave-one-glacier-out cross-validation (N = 255). Left: error distribution of the computed mass-balance if determined by the interpolated \(t^*\). Right: error distribution of the mass-balance if determined by interpolation of \(\mu^*\). From Maussion et al., 2019.

If you come up with a great idea to make the calibration procedure better, please reach out! We have some ideas, but haven’t come up to try it yet. There are so many things we want to try out with OGGM!

Prepared reference t* lists

As of OGGM v1.4, we offer various t* lists for download and use in your workflow. You will find them here:


The easiest is to let OGGM download them for you with, for example:

# Reference data for the recalibration of the mass-balance
ref_tstars_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5'

See also: local_t_star documentation.

Bonus: additional MB data for OGGM

This is a rare, but possible use case: you might have mass-balance data that is not yet available in the WGMS database, or reanalyzed data that you would like to use in place of the default WGMS. Doing this is quite simple: let’s assume we have some secret data for a third glacier, here the Oberer Grindelwald Glacier in the Swiss Alps (this is an example of course):

gdirs = workflow.init_glacier_directories(['RGI60-11.00787', 'RGI60-11.00897', 'RGI60-11.01270'], from_prepro_level=3, reset=True, force=True)
# Let's make sure that OGG is not a reference glacier: this line would throw an error:
# gdirs[2].get_ref_mb_data()

First, let’s start by adding some data to this glacier. We are just creating some fake data and giving it to the glacier directory object.

# This dataframe just needs some index and an 'ANN'
df = pd.DataFrame(index=range(1980, 2010))
# Get random but repeatable results
df['ANNUAL_BALANCE'] = np.random.randn(len(df.index)) * 800 - 200

# Give it to the object
gdir = gdirs[2]  # get OGG

Second, we have to tell OGGM that this glacier is a reference glacier. We do this by adding it’s ID to the list:

# Make sure the default list is generated
# Append our glacier
# Make sure all three glaciers are now reference glaciers
assert len(utils.get_ref_mb_glaciers(gdirs)) == 3

Third, re-run the calibration!

cfg.PARAMS['run_mb_calibration'] = True
workflow.execute_entity_task(tasks.local_t_star, gdirs);
workflow.execute_entity_task(tasks.mu_star_calibration, gdirs);

Note: what just happened here? We asked OGGM to recalibrate its mass-balance model. In essence, it is as simple as computing a list of so-called “reference \(t^*\)” from which the model can interpolate from. This list is then written as CSV file in the working directory:

df_tstar = pd.read_csv(os.path.join(cfg.PATHS['working_dir'], 'ref_tstars.csv'), index_col=0)

This list has to be generated once, and has to be copied in each working directory you use if you want OGGM to take it into account. In the absence of this file in the working directory, OGGM will use its default pre-calibrated file, located in the online sample data repository.

After copying the file in a fresh working directory, the task compute_ref_t_stars should not be called. Only local_t_star and mu_star_calibration matter.

Let’s see how well our new glacier’s mass-balance can be reproduced:

ref_df = gdir.get_ref_mb_data()
mbmod = massbalance.PastMassBalance(gdir)
ref_df['OGGM (fake calib)'] = mbmod.get_specific_mb(fls=gdir.read_pickle('inversion_flowlines'), 
ref_df[['ANNUAL_BALANCE', 'OGGM (fake calib)']].mean()

A bias of zero! Works perfectly, right? ;-)

ref_df[['ANNUAL_BALANCE', 'OGGM (fake calib)']].plot();

Take home points

  • We illustrated with a simplified example how to do a simple cross-validation of the mass-balance calibration of OGGM. In reality, the \(t^*\) are interpolated, not assigned like this. The results of the full cross-validation are therefore different for Kesselwandferner.

  • It is possible to play around with OGGM mass-balance data for sensitivity experiments.

  • It is also possible to add custom mass-balance data to OGGM glaciers and re-calibrate the model based on this data

  • The mass-balance model of OGGM can lead to substantial biases at the local scale. We show that on average the bias is low, but individual glaciers can have substantial biases.

  • Future improvements will have to deal with this bias, e.g. by making use of the increasingly available geodetic mass-balance measurements.

What’s next?