# A look into the new mass balance calibration in OGGM v1.6#

The default mass-balance (MB) model of OGGM is a very standard temperature index melt model.

In versions before 1.6, OGGM had a complex calibration procedure which originated from the times where we had only observations from a few hundred glaciers. We used them to calibrate the model and then a so-called tstar (famous in only very niche circles) which was interpolated to glaciers without observations (see the original publication). This method was very powerful but, as new observational datasets emerged, we can now calibrate on a glacier-per-glacier basis. With the new era of geodetic observations, OGGM uses per default the average geodetic observations from Jan 2000–Jan 2020 of Hugonnet al. 2021, that are now available for almost every glacier world-wide.

Pre-processed directories from OGGM (from the Bremen server) have been calibrated for you, based on a specific climate dataset (W5E5) and our own dedicated calibration strategy. But, what if you want to use another climate dataset? Or another reference dataset?

In this notebook, we will:

• introduce the new flexible calibration method available to all users since OGGM v1.6: mb_calibration_from_scalar_mb

• show how users can apply it (and extend it) for their own workflow

• explain how we used it to calibrate the OGGM pre-processed directories

## Set-up#

import matplotlib.pyplot as plt
import matplotlib
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
from oggm.core.massbalance import mb_calibration_from_scalar_mb, mb_calibration_from_geodetic_mb, mb_calibration_from_wgms_mb

cfg.initialize(logging_level='WARNING')
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-calib-mb', reset=True)
cfg.PARAMS['border'] = 10

2023-03-10 19:21:23: oggm.cfg: Reading default parameters from the OGGM params.cfg configuration file.
2023-03-10 19:21:23: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-10 19:21:23: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-03-10 19:21:23: oggm.cfg: PARAMS['border'] changed from 80 to 10.


We start from two well known glaciers in the Austrian Alps, Kesselwandferner and Hintereisferner. But you can also choose any other other glacier, e.g. from this list.

# we start from preprocessing level 3
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.1/elev_bands/W5E5/'
gdirs = workflow.init_glacier_directories(['RGI60-11.00787',  'RGI60-11.00897'], from_prepro_level=3, prepro_base_url=base_url)

2023-03-10 19:21:24: oggm.workflow: init_glacier_directories from prepro level 3 on 2 glaciers.
2023-03-10 19:21:24: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 2 glaciers


The two glaciers are neighbors but have very different geometries:

f, ax = plt.subplots(1,1,figsize=(6, 6))
ax.set_title(gdirs[0].rgi_id + ' & ' + gdirs[1].rgi_id);


## Calibrated mass balance parameters in OGGM v1.6 pre-processed directories#

We just downloaded the data. Let’s have a look at the climate data we used and the calibrated parameters:

# Pick each glacier
gdir_kwf = gdirs[0]
gdir_hef = gdirs[1]

gdir_hef.get_climate_info()

{'baseline_climate_source': 'GSWP3_W5E5',
'baseline_yr_0': 1901,
'baseline_yr_1': 2019,
'baseline_climate_ref_hgt': 2252.0,
'baseline_climate_ref_pix_lon': 10.75,
'baseline_climate_ref_pix_lat': 46.75}


This is called the “baseline climate” in OGGM and is necessary to calibrate the model against observations. Ideally, the baseline climate should be real observations as perfect as possible, but in reality this is not the case. Often, gridded climate datasets have biases - we need to take this into accound during our calibration. Let’s have a look at the mass balance parameters for both glaciers:

gdir_hef.read_json('mb_calib')

{'rgi_id': 'RGI60-11.00897',
'bias': 0,
'melt_f': 5.0,
'prcp_fac': 3.532479282420574,
'temp_bias': 1.758147556769245,
'reference_mb': -1100.3,
'reference_mb_err': 171.8,
'reference_period': '2000-01-01_2020-01-01',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}

gdir_kwf.read_json('mb_calib')

{'rgi_id': 'RGI60-11.00787',
'bias': 0,
'melt_f': 5.0,
'prcp_fac': 2.96302677905353,
'temp_bias': 1.758147556769245,
'reference_mb': -514.8000000000001,
'reference_mb_err': 136.3,
'reference_period': '2000-01-01_2020-01-01',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}


We will explain later what all these mean. Lets focus on these for now: melt_f, prcp_fac, and temp_bias which have been calibrated with reference_mb over the reference period reference_period.

Per default the Hugonnet et al. (2021) average geodetic observation is used over the entire time period Jan 2000 to Jan 2020 to calibrate the MB model parameter(s) for every single glacier.

For our two neighboring example glaciers, that share the same climate gridpoint, the same pre-calibrated temp_bias and the same melt_f is used. The prcp_fac is slighly different, hence in that example, changing the prcp_fac within the narrow range of $$[0.8,1.2]$$ was sufficient to match the MB model to the observations.

Note that the two glaciers are in the same climate (from the forcing data) but are very different in size, orientation and geometry. So, at least one of the MB model parameters is different (here the prcp_fac) and also the observed MB values vary:

for gdir in gdirs:
mbmod = massbalance.MonthlyTIModel(gdir)
year=np.arange(2000,2020,1)).mean()
print(gdir.rgi_id, f': average MB 2000-2020 = {mean_mb:.1f} kg m-2, ',

RGI60-11.00787 : average MB 2000-2020 = -514.8 kg m-2,  prcp_fac: 2.96
RGI60-11.00897 : average MB 2000-2020 = -1100.3 kg m-2,  prcp_fac: 3.53


Kesselwandferner has a less negative MB than its neighbor, probably because it is smaller in size and spans less altitude. To simplify things, we will now focus just on one glacier (Hintereisferner), but you can repeat it of course with other glaciers.

How are these parameters calibrated by OGGM? How can I calibrate them for my use case?
These are the questions we want to address in the rest of the notebook. We will change the three MB model parameters (melt_f, temp_bias and prcp_fac) and see their influence on the calibration, to allow users to make an informed decision on how to calibrate them for their purpose. Then we will go back to the calibration method used by OGGM for the preprocessed directories.

Note alst that there are some global MB parameters (mb_global_params), which we assume to be the same globally for every glacier:

gdir_hef.read_json('mb_calib')['mb_global_params']

{'temp_default_gradient': -0.0065,
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0}


These global MB parameter values were found to represent best the in-situ observations during a cross-validation for glaciers with additional observations. Of course, they could also be changed for different glaciers, but this would need even more data to justify the differences! The influence of using a different temperature lapse rate gradient to the default -0.0065 K/km are analysed in Schuster et al. (2023, in review).

First, let’s check if the OGGM calibration worked at all

• we will get first the average modelled MB:

h, w = gdir_hef.get_inversion_flowline_hw()
mb_geod = massbalance.MonthlyTIModel(gdir_hef)
# Note, if you change cfg.PARAMS['geodetic_mb_period'], you need to change the years here as well!
mbdf= pd.DataFrame(index = np.arange(2000,2020,1))
mbdf['mod_mb'] = mb_geod.get_specific_mb(h, w, year=mbdf.index)
mbdf.mean()

mod_mb   -1100.3
dtype: float64

• then get the observed geodetic MB that we calibrated our mass-balance to:

print('reference MB: ' + str(gdir_hef.read_json('mb_calib')['reference_mb'])+ ' kg m-2')
print('reference MB uncertainties: '+ str(gdir_hef.read_json('mb_calib')['reference_mb_err'])+ ' kg m-2')
print('reference MB time period: ' + gdir_hef.read_json('mb_calib')['reference_period'])

reference MB: -1100.3 kg m-2
reference MB uncertainties: 171.8 kg m-2
reference MB time period: 2000-01-01_2020-01-01

• We calibrated on the entire 20-yr time period, because then the observational uncertainties are smallest.

# this tests if the two parameters are very similar:


Perfect! Our MB model reproduces the average observed MB, so the calibrated worked!

## The new flexible calibration scheme in OGGM#

We have now used the option that is applied in the preprocessed directories. However, there are many other calibration options. We will show you some, but you could also invent your own for your purposes.

Calibrate melt_f, fixed prcp_fac and temp_bìas

For example, let’s calibrate the melt_f to match the same observations as before (i.e., average geodetic observation) and fix the other parameters. We will use the more general function mb_calibration_from_scalar_mb which is most flexible. However, it requires you to give manually the observed MB on which you want to calibrate:

# if you use the default calibration data from Hugonnet et al., 2021,
# we can get the geodetic data from any glacier from here:
ref_mb_df = utils.get_geodetic_mb_dataframe().loc[gdir_hef.rgi_id]
ref_mb_df = ref_mb_df.loc[ref_mb_df['period'] == cfg.PARAMS['geodetic_mb_period']]
# dmdtda: in meters water-equivalent per year -> we convert to kg m-2 yr-1
ref_mb = float(ref_mb_df['dmdtda']) * 1000
ref_mb

-1100.3

# Let's calibrate the melt_f (this is the default option in mb_calibration_from_scalar_mb)
# overwrite_gdir has to be set to True,
# because we want to overwrite the old calibration
mb_calibration_from_scalar_mb(gdir_hef,
ref_mb = ref_mb,
ref_period=cfg.PARAMS['geodetic_mb_period'],
overwrite_gdir=True);

mb_melt_f = massbalance.MonthlyTIModel(gdir_hef)
mbdf['mod_mb_melt_f'] = mb_melt_f.get_specific_mb(h, w, year=mbdf.index)


Let’s read the new calibrated MB model parameter for the glacier:

gdir_hef.read_json('mb_calib')

{'rgi_id': 'RGI60-11.00897',
'bias': 0,
'melt_f': 9.06088881579363,
'prcp_fac': 3.3571363162713674,
'temp_bias': 0,
'reference_mb': -1100.3,
'reference_mb_err': None,
'reference_period': '2000-01-01_2020-01-01',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}


Only melt_f is used to calibrate the MB model. prcp_fac was chosen depending on the average winter precipitation (explained below) and temp_bias is 0.

Calibrate temp_bias, fixed prcp_fac and melt_f

We will now instead fix the melt_f and the prcp_fac, and use the temp_bias as free variable for calibration. For that we overwrite the default value of calibrate_param1:

# Let's calibrate on the temp_bias instead
# overwrite_gdir has to be set to True,
# because we want to overwrite the old calibration
mb_calibration_from_scalar_mb(gdir_hef,
ref_mb = ref_mb,
ref_period=cfg.PARAMS['geodetic_mb_period'],
calibrate_param1='temp_bias',
overwrite_gdir=True)

mb_temp_b = massbalance.MonthlyTIModel(gdir_hef)
mbdf['mod_mb_temp_b'] = mb_temp_b.get_specific_mb(h, w, year=mbdf.index)

mb_params = gdir_hef.read_json('mb_calib')
mb_params['temp_bias'], mb_params['prcp_fac'], mb_params['melt_f']

(1.6621867397806114, 3.3571363162713674, 5.0)


Here we used the median melt_f (i.e., 5 kg m-2 day-1 K-1) and the same prcp_fac as in the previous option, but changed temp_bias until the ref_mb is matched. As the melt_f is lower than in the previous calibration option, we need to have a positive temp_bias to get to the same average MB (i.e., ref_mb).

Calibrate prcp_fac, fixed temp_bias and melt_f

Similarly, we can also just calibrate on the precipitation factor (prcp_fac).

# Let's calibrate on the prcp_fac instead
# overwrite_gdir has to be set to True,
# because we want to overwrite the old calibration
mb_calibration_from_scalar_mb(gdir_hef,
ref_mb = ref_mb,
ref_period=cfg.PARAMS['geodetic_mb_period'],
calibrate_param1='prcp_fac',
overwrite_gdir=True)

mb_prcp_fac = massbalance.MonthlyTIModel(gdir_hef)
mbdf['mod_mb_prcp_fac'] = mb_prcp_fac.get_specific_mb(h, w, year=mbdf.index)

mb_params['temp_bias'], mb_params['prcp_fac'], mb_params['melt_f']

(0, 1.0160354839620402, 5.0)


We chose two glaciers that actually have also in-situ observations available. We will get the in-situ observations like that:

mbdf['in_situ_mb'] = gdir_hef.get_ref_mb_data().loc[2000:2019]['ANNUAL_BALANCE']
# if we want to compare the average to the average geodetic MB, we need to have all years available
assert len(mbdf['in_situ_mb']) == 20


Let’s plot how well we match the interannual observations for the different calibration options:

plt.plot(mbdf['in_situ_mb'], label='in-situ observations\n'+f'average 2000-2020: {mbdf.in_situ_mb.mean():.1f} '+ r'kg m$^{-2}$', color='k', lw=3)
plt.plot(mbdf['mod_mb'],
label='modelled mass-balance via preprocessed directories\n'+f'average 2000-2020: {mbdf.mod_mb.mean():.1f} ' + r'kg m$^{-2}$')
plt.plot(mbdf['mod_mb_melt_f'],
label='modelled mass-balance via calibrating melt_f\n'+f'average 2000-2020: {mbdf.mod_mb.mean():.1f} ' + r'kg m$^{-2}$')
plt.plot(mbdf['mod_mb_temp_b'],
label='modelled mass-balance via calibrating temp_bias\n'+f'average 2000-2020: {mbdf.mod_mb_temp_b.mean():.1f} ' + r'kg m$^{-2}$')
plt.plot(mbdf['mod_mb_prcp_fac'],
label='modelled mass-balance via calibrating prcp_fac\n'+f'average 2000-2020: {mbdf.mod_mb_prcp_fac.mean():.1f} ' + r'kg m$^{-2}$')
plt.xticks(np.arange(2000,2020,2))
plt.legend(bbox_to_anchor=(1,1))
plt.ylabel(r'specific mass-balance (kg m$^{-2}$)')
plt.xlabel('Year');


For this glacier, over the same time period (here 2000-2020), the average MB is quite similar between the geodetic and in-situ observation and could even be explained by the fact that the in-situ observations are in hydrological years (starting in October of the previous year) and the geodetic observations are in calendar years. If you repeat the analysis for other glaciers with in-situ observations over that time period, you might see larger discrepancies.

You can also see, that there are differences in the annual MB between the calibration options and the in-situ observations. We will analyse these differences more systematically later!

## Calibrate on direct glaciological in-situ MB observations#

Ok, but actually you might trust the in-situ observations much more than the geodetic observation. So if you are only interested in glaciers with these in-situ observations, you can also use the in-situ observations for calibration. This might allow you also to calibrate over a longer time period and to use additional informations (such as the interannual MB variability, seasonal MB, …). However, bare in mind that there are often gaps in the in-situ MB time series and make sure that you calibrate to the same modelled time period!

Attention: For the Hintereisferner glacier with an almost 70-year long time series, the assumption that we make, i.e., that the area does not change over the time period, gets more problematic. So think twice before repeating this at home!

mb_calibration_from_wgms_mb(gdir_hef, overwrite_gdir=True)

{'rgi_id': 'RGI60-11.00897',
'bias': 0,
'melt_f': 11.27883598160852,
'prcp_fac': 3.3571363162713674,
'temp_bias': 0,
'reference_mb': -648.223880597015,
'reference_mb_err': None,
'reference_period': 'custom',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}


mb_calibration_from_wgms_mb is a short-cut function that uses internally mb_calibration_from_scalar_mb, but applies directly the average annual in-situ MB observations from the WGMS. Per default, the melt_f is calibrated, and prcp_fac and temp_bias are fixed. You can get the observations like that:

mbdf_in_situ = gdir_hef.get_ref_mb_data()
mbdf_in_situ['ref_mb'] = mbdf_in_situ['ANNUAL_BALANCE']
# check that every year between the beginning and the end has MB observations
assert len(mbdf_in_situ.index) == (mbdf_in_situ.index[-1] - mbdf_in_situ.index[0] + 1)
ref_mb = mbdf_in_situ.ref_mb.mean()

mb_in_situ_obs = massbalance.MonthlyTIModel(gdir_hef)
mbdf_in_situ['mod_mb'] = mb_in_situ_obs.get_specific_mb(h, w, year=mbdf_in_situ.index)
plt.plot(mbdf_in_situ['ref_mb'],
label='in-situ observations\n'+f'average: {ref_mb:.1f} '+ r'kg m$^{-2}$', color='grey', lw=3)
plt.plot(mbdf_in_situ['mod_mb'],
label=('modelled mass-balance\nvia calibrating melt_f\n'
+f'average: {mbdf_in_situ.mod_mb.mean():.1f} ' + r'kg m$^{-2}$'))

plt.legend()
plt.ylabel(r'specific mass-balance (kg m$^{-2}$)')
plt.xlabel('Year');


The average MB over the entire time period is matched as we calibrate to it. However, the interannual mass-balance variability is different. How well the modelled annual mass-balance time series matches the observations can be for example assessed by looking at the correlation between modelled and observed annual MB

mbdf_in_situ[['ref_mb','mod_mb']].corr().values[0][1]

0.806883065989496


or by looking into the standard deviation of interannual MB variability, where we see a much larger variability for the modelled MB:

mbdf_in_situ[['ref_mb','mod_mb']].std()

ref_mb     624.538759
mod_mb    1276.246655
dtype: float64


### Dealing with errors and including your own mass-balance observations#

Let’s look at another glacier, now in Central Asia, for example the Golubin glacier:

gdir_gol = workflow.init_glacier_directories(['RGI60-13.11609'], from_prepro_level=3, prepro_base_url=base_url)[0]

f, ax = plt.subplots(1,1,figsize=(6, 6))
plt.tight_layout()

2023-03-10 19:21:29: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2023-03-10 19:21:29: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers


Let’s calibrate that glacier by only changing the melt_f, default in mb_calibration_from_scalar_mb:

# first need to get the MB observations again
ref_mb_df = utils.get_geodetic_mb_dataframe().loc[gdir_gol.rgi_id]
ref_mb_df = ref_mb_df.loc[ref_mb_df['period'] == cfg.PARAMS['geodetic_mb_period']]
# dmdtda: in meters water-equivalent per year -> we convert to kg m-2 yr-1
ref_mb = float(ref_mb_df['dmdtda']) * 1000
ref_mb

-255.9

# Note that this will create an error!
mb_calibration_from_scalar_mb(gdir_gol, ref_mb=ref_mb, ref_period = cfg.PARAMS['geodetic_mb_period'])

2023-03-10 19:21:32: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-13.11609: RGI60-13.11609: ref mb not matched. Try to set calibrate_param2.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/core/massbalance.py:1778, in mb_calibration_from_scalar_mb(gdir, ref_mb, ref_mb_err, ref_period, ref_mb_years, write_to_gdir, overwrite_gdir, calibrate_param1, calibrate_param2, calibrate_param3, melt_f, melt_f_min, melt_f_max, prcp_fac, prcp_fac_min, prcp_fac_max, temp_bias, temp_bias_min, temp_bias_max, mb_model_class)
1777 try:
-> 1778     optim_param1 = optimize.brentq(to_minimize,
1779                                    min_range, max_range,
1780                                    args=(calibrate_param1,)
1781                                    )
1782 except ValueError:

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/scipy/optimize/_zeros_py.py:783, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp)
782     raise ValueError("rtol too small (%g < %g)" % (rtol, _rtol))
--> 783 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
784 return results_c(full_output, r)

ValueError: f(a) and f(b) must have different signs

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
Cell In[29], line 2
1 # Note that this will create an error!
----> 2 mb_calibration_from_scalar_mb(gdir_gol, ref_mb=ref_mb, ref_period = cfg.PARAMS['geodetic_mb_period'])

490 ex_t = time.time()
--> 491 out = task_func(gdir, **kwargs)
492 ex_t = time.time() - ex_t

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/core/massbalance.py:1784, in mb_calibration_from_scalar_mb(gdir, ref_mb, ref_mb_err, ref_period, ref_mb_years, write_to_gdir, overwrite_gdir, calibrate_param1, calibrate_param2, calibrate_param3, melt_f, melt_f_min, melt_f_max, prcp_fac, prcp_fac_min, prcp_fac_max, temp_bias, temp_bias_min, temp_bias_max, mb_model_class)
1782 except ValueError:
1783     if not calibrate_param2:
-> 1784         raise RuntimeError(f'{gdir.rgi_id}: ref mb not matched. '
1785                            f'Try to set calibrate_param2.')
1787     # Check which direction we need to go
1788     diff_1 = to_minimize(min_range, calibrate_param1)

RuntimeError: RGI60-13.11609: ref mb not matched. Try to set calibrate_param2.


We got a RuntimeError that says that the ref_mb is not matched and that we should set calibrate_param2. What happened?

Well, no melt_f parameter could be found in the given ranges to reproduce the given calibration data (ref_mb).

# What are the current minimum and maximum ranges of the melt factor (unit: kg m-2 day-1 K-1)
cfg.PARAMS['melt_f_min'], cfg.PARAMS['melt_f_max']

(1.5, 17.0)


If we believe that our observations are true, maybe the climate or the processes represented by the MB model are erroneous.

What can we do to still match the ref_mb? We can either change the melt_f ranges by setting other values to the parameters above or we can allow that another MB model parameter is changed (i.e., calibrate_param2). This is basically very similar to the three-step-calibration first introduced in Huss & Hock 2015, but here you can choose your parameter ranges and parameter order yourself.

To reduce these MB model calibration errors, you can first change the melt_f, fix it at the lower or upper limit, and then change the temp_bias.

# Allowing another parameter to change is done by defining calibrate_param2
mb_calibration_from_scalar_mb(gdir_gol, ref_mb=ref_mb,
write_to_gdir=False,
ref_period =cfg.PARAMS['geodetic_mb_period'],
calibrate_param1='melt_f',
calibrate_param2='temp_bias')

{'rgi_id': 'RGI60-13.11609',
'bias': 0,
'melt_f': 1.5,
'prcp_fac': 3.7686569086250965,
'temp_bias': -0.5558387442516765,
'reference_mb': -255.9,
'reference_mb_err': None,
'reference_period': '2000-01-01_2020-01-01',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}


The melt_f is at the minimum value and temperature was corrected to more negative values to allow for the relatively weak negative MB.

Use your own or fake MB observations for calibration

In the same way, you can also use your own MB observations or fake data to calibrate the MB model. We use here as an example, an unrealistically high positive MB for the Hintereisferner over a 10-year time period. To allow the calibration to happen, we need again to set calibrate_param2! Here we will use the prcp_fac as second parameter for a change:

ref_period = '2000-01-01_2010-01-01'
ref_mb = 2000 # Let's use an unrealistically positive  mass-balance
mb_calibration_from_scalar_mb(gdir_hef, ref_mb=ref_mb,
ref_period=ref_period,
overwrite_gdir=True,
calibrate_param2='prcp_fac');

mb_new = massbalance.MonthlyTIModel(gdir_hef)
# ok, we actually matched the new ref_mb over the 10-yr period
np.testing.assert_allclose(ref_mb,
mb_new.get_specific_mb(h, w, year=np.arange(2000,2010,1)).mean())
# Let's look at the calibrated parameters

{'rgi_id': 'RGI60-11.00897',
'bias': 0,
'melt_f': 1.5,
'prcp_fac': 4.217210907878396,
'temp_bias': 0,
'reference_mb': 2000,
'reference_mb_err': None,
'reference_period': '2000-01-01_2010-01-01',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}


Ok, in that case, the climate was too warm to allow for such a positive MB. Even the smallest melt_f did not allow for such a positive MB. Therefore, the prcp_fac was increased to allow that more (solid) precipitation can occur.

However, also the prcp_fac has a limited range:

cfg.PARAMS['prcp_fac_min'], cfg.PARAMS['prcp_fac_max']

(0.1, 10.0)


So, if we increase the ref_mb even further, we might even need a third free parameter (i.e., we need to set calibrate_param3). Let’s try it out, by using now the same order as in Huss & Hock (2015) (but different allowed parameter changes that also depend on the chosen climate…):

That means, we match the ref_mb by

• first trying to adapt prcp_fac,

• then melt_f,

• and finally temp_bias.

ref_mb = 6000 # if you replace it with 8000, no combination of parameters can be found for that glacier
mb_calibration_from_scalar_mb(gdir_hef,ref_mb=ref_mb,
ref_period=ref_period,
calibrate_param1='prcp_fac',
calibrate_param2='melt_f',
calibrate_param3='temp_bias',
overwrite_gdir=True)

{'rgi_id': 'RGI60-11.00897',
'bias': 0,
'melt_f': 1.5,
'prcp_fac': 10.0,
'temp_bias': -0.8132605320838399,
'reference_mb': 6000,
'reference_mb_err': None,
'reference_period': '2000-01-01_2010-01-01',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}


In that case, the minimum melt_f and the maximum prcp_fac are applied. The temp_bias is set to a negative value as this results in more solid precipitation. If you increased the ref_mb even further, the method will not find any combination as the temp_bias also has a limited range. If you really want to match the observation, then you would need to change the parameter ranges.

cfg.PARAMS['temp_bias_min'], cfg.PARAMS['temp_bias_max']

(-15.0, 15.0)

### check out the docstring for more information about the options
### by uncommenting the line below:
# mb_calibration_from_scalar_mb?


## Overparameteristion or the magic choice of the best calibration option:#

We found already some combinations that equally well match the average MB over a given time period. As we only use only one observation per glacier (i.e., per default the average geodetic MB from 2000-2020), but have up to three free MB model parameters, the MB model is overparameterised. That means, there are in theory an infinite amount of calibration options possible that equally well match the one obervation. Let’s look a bit more systematically into that:

We will use a range of different prcp_fac and then calibrate the melt_f accordingly to always match to the default average MB (ref_mb) over the reference period (ref_period).

# let's go back to the Hintereisferner glacier and get again the MB observations:
ref_mb_df = utils.get_geodetic_mb_dataframe().loc[gdir_hef.rgi_id]
ref_mb_df = ref_mb_df.loc[ref_mb_df['period'] == cfg.PARAMS['geodetic_mb_period']]
ref_mb = float(ref_mb_df['dmdtda']) * 1000

# calibrate the melt_f and annual MB
pd_prcp_fac_sens = pd.DataFrame(index=np.arange(0.1,5.0,0.3))
spec_mb_prcp_fac_sens_dict = {}
for prcp_fac in pd_prcp_fac_sens.index:
calib_param = mb_calibration_from_scalar_mb(gdir_hef, ref_mb=ref_mb, ref_period=cfg.PARAMS['geodetic_mb_period'],
prcp_fac=prcp_fac, overwrite_gdir=True)
pd_prcp_fac_sens.loc[prcp_fac, 'melt_f'] = calib_param['melt_f']
mb_prcp_fac_sens = massbalance.MonthlyTIModel(gdir_hef)
# ok, we actually matched the new ref_mb
annual_mb = mb_prcp_fac_sens.get_specific_mb(h, w, year=np.arange(2000,2020,1))
spec_mb_prcp_fac_sens_dict[prcp_fac] = annual_mb
# let's also check how well the interannual MB variability is matched
# by comparing the standard deviation of the annual MB to the observed one
std_quot_annual_mb = annual_mb.std()/mbdf_in_situ.loc[2000:2019, 'ANNUAL_BALANCE'].std()
pd_prcp_fac_sens.loc[prcp_fac, 'std_quot_annual_mb'] = std_quot_annual_mb

# let's get some nice colors visualising more precipitation
colors_prcp_fac = plt.get_cmap('viridis_r').colors[10::10]
plt.figure()
for j,prcp_fac in enumerate(pd_prcp_fac_sens.index):
plt.plot(prcp_fac, pd_prcp_fac_sens.loc[prcp_fac, 'melt_f'],
'o', color=colors_prcp_fac[j])
plt.ylabel(r'melt_f (kg m$^{-2}$ day$^{-1}$ K$^{-1}$)')
plt.xlabel('prcp_fac');


All these combinations match the average ref_mb equally. The larger the chosen prcp_fac, the larger needs to be the calibrated melt_f when matching to the same average MB.

What is the influence on the chosen parameter combination on other estimates than the average MB?

plt.figure()
for j,prcp_fac in enumerate(pd_prcp_fac_sens.index):
plt.plot(np.arange(2000,2020,1),
spec_mb_prcp_fac_sens_dict[prcp_fac], '-',
color=colors_prcp_fac[j], label=f'{prcp_fac:.1f}')
plt.plot(mbdf_in_situ.loc[2000:2019].index,
mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'],
color='grey', lw=3,
label='observed in-situ')
plt.ylabel(r'Annual mass-balance (kg m$^{-2}$)')
plt.xlabel('Year')
plt.legend(title='prcp_fac:', bbox_to_anchor=(1,1), fontsize=9)
plt.xticks(np.arange(2000,2020,2));
plt.title(gdir_hef.rgi_id);


The larger the prcp_fac and the melt_f, the larger is the interannual MB variability. For glaciers with in-situ observations, we can find a combination of prcp_fac and melt_f that has a similar interannnual MB variability than the observations. For example, you can choose a MB model parameter combination where the standard deviation quotient of the annual the modelled and observed MB is near to 1:

condi = (pd_prcp_fac_sens['std_quot_annual_mb'] < 1.1) & (pd_prcp_fac_sens['std_quot_annual_mb'] > 0.9)
pd_prcp_fac_sens[condi]

melt_f std_quot_annual_mb
1.6 6.012949 0.952705
1.9 6.533331 1.052461

For the Hintereisferner glacier, the best MB model combination to match the interannual MB variability would be for a prcp_fac between 1.6 and 1.9 and a melt_f between 6 and 6.5 (more in Schuster et al., 2023, in review).

We can also fix the prcp_fac, and change the temp_bias and melt_f instead:

melt_f_dict_tb = {}
spec_mb_temp_b_sens_dict = {}

for temp_bias in np.arange(-5,5.0,0.5):
# for too negative temp_bias, no melt_f is found that matches the observations. We would need to
# change the prcp_fac , but here we will just look
# at those combinations where calibration works with a fixed prcp_fac.
try:
calib_param = mb_calibration_from_scalar_mb(gdir_hef, ref_mb=ref_mb, ref_period=cfg.PARAMS['geodetic_mb_period'],
temp_bias=temp_bias, overwrite_gdir=True)
melt_f_dict_tb[temp_bias] = calib_param['melt_f']
mb_temp_b_sens = massbalance.MonthlyTIModel(gdir_hef)
# ok, we actually matched the new ref_mb
spec_mb_temp_b_sens_dict[temp_bias] = mb_temp_b_sens.get_specific_mb(h, w, year=np.arange(2000,2020,1))
except RuntimeError: #, 'RGI60-11.00897: ref mb not matched. Try to set calibrate_param2'
pass

2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:35: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.

# let's get a nice colormap centered at temp_bias=0
norm = matplotlib.colors.Normalize(vmin=-5, vmax=5.01)
colors_temp_bias = plt.get_cmap('coolwarm')

fig, axs = plt.subplots(1,2,figsize=(14,6))
for j,temp_bias in enumerate(melt_f_dict_tb.keys()):
axs[0].plot(temp_bias, melt_f_dict_tb[temp_bias], 'o',
color=colors_temp_bias(norm(temp_bias)))
axs[0].set_ylabel(r'melt_f (kg m$^{-2}$ day$^{-1}$ K$^{-1}$)')
axs[0].set_xlabel('temp_bias (°C)');

for temp_bias in melt_f_dict_tb.keys():
axs[1].plot(np.arange(2000,2020,1),
spec_mb_temp_b_sens_dict[temp_bias], '-',
color=colors_temp_bias(norm(temp_bias)),
label=temp_bias)
axs[1].plot(mbdf_in_situ.loc[2000:2019].index, mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'], color='grey', lw=3,
label='observed in-situ')
axs[1].set_ylabel(r'Annual mass-balance (kg m$^{-2}$)')
axs[1].set_xlabel('Year')
axs[1].legend(title='temp_bias:', bbox_to_anchor=(1,1))
axs[1].set_xticks(np.arange(2000,2020,2));
axs[1].set_title(gdir_hef.rgi_id);
plt.tight_layout()


A lower melt_f is needed if a positive temp_bias is applied… and with that the interannual MB variability gets smaller.

Or we change prcp_fac and temp_bias, and fix the melt_f:

pf_temp_b_dict_tb = {}
spec_mb_pf_temp_b_sens_dict = {}

for temp_bias in np.arange(-5,5.0,0.5):
# for too negative temp_bias, no prcp_fac is found that matches the observations. We would need to
# change the melt_f , but here we will just look
# at those combinations where calibration works with a fixed melt_f.
try:
calib_param = mb_calibration_from_scalar_mb(gdir_hef, calibrate_param1='prcp_fac',
ref_mb=ref_mb, ref_period=cfg.PARAMS['geodetic_mb_period'],
temp_bias=temp_bias, overwrite_gdir=True)
pf_temp_b_dict_tb[temp_bias] = calib_param['prcp_fac']
mb_pf_temp_b_sens = massbalance.MonthlyTIModel(gdir_hef)
# ok, we actually matched the new ref_mb
spec_mb_pf_temp_b_sens_dict[temp_bias] = mb_pf_temp_b_sens.get_specific_mb(h, w, year=np.arange(2000,2020,1))
except RuntimeError: #, 'RGI60-11.00897: ref mb not matched. Try to set calibrate_param2'
pass

2023-03-10 19:21:36: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:36: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:36: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:37: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:37: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:37: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:37: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:37: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.
2023-03-10 19:21:37: oggm.core.massbalance: RuntimeError occurred during task mb_calibration_from_scalar_mb on RGI60-11.00897: RGI60-11.00897: ref mb not matched. Try to set calibrate_param2.

# let's get another centered colormap at temp_bias=0
norm = matplotlib.colors.Normalize(vmin=-5, vmax=5.01)
colors_temp_bias = plt.get_cmap('Spectral_r')

fig, axs = plt.subplots(1,2,figsize=(14,6))
for j,temp_bias in enumerate(pf_temp_b_dict_tb.keys()):
axs[0].plot(temp_bias, pf_temp_b_dict_tb[temp_bias], 'o',
color=colors_temp_bias(norm(temp_bias)))
axs[0].set_ylabel(r'prcp_fac')
axs[0].set_xlabel('temp_bias (°C)');

for temp_bias in pf_temp_b_dict_tb.keys():
axs[1].plot(np.arange(2000,2020,1),
spec_mb_pf_temp_b_sens_dict[temp_bias], '-',
color=colors_temp_bias(norm(temp_bias)),
label=temp_bias)
axs[1].plot(mbdf_in_situ.loc[2000:2019].index, mbdf_in_situ.loc[2000:2019]['ANNUAL_BALANCE'], color='grey', lw=3,
label='observed in-situ')
axs[1].set_ylabel(r'Annual mass-balance (kg m$^{-2}$)')
axs[1].set_xlabel('Year')
axs[1].legend(title='temp_bias:', bbox_to_anchor=(1,1))
axs[1].set_xticks(np.arange(2000,2020,2));
axs[1].set_title(gdir_hef.rgi_id);
plt.tight_layout()


A larger temp_bias needs a larger prcp_fac … and with that the interannual MB variability gets larger.

The parameter combination choice or calibration option also influences the modelled seasonal MB or elevation-dependent MB over the calibration period. Additionally, volume and runoff are influenced by the model parameter choice. If you are further interested in that, you can have a look into Schuster et al (2023, in review), which analyses the “Glacier projections sensitivity to temperature-index model choices and calibration strategies”.

Using prior knowledge on the MB model parameters by a Bayesian calibration scheme is a good option to constrain the parameter ranges. Like that you can also directly include the MB observation uncertainties and quantify the parameter uncertainties. You can read about that in Rounce et al., 2020!

We are working on including a Bayesian calibration framework into OGGM by using the calibration procedure explained below as prior estimates.

## Back to where we started: how did we calibrate the default pre-processed directories in OGGM 1.6?#

There are no easy solutions to the problem of overparameterisation. Here we explain the steps taken by OGGM to calibrate the parameters globally for our users who don’t want to care about calibration.

In the preprocessed directories MB model calibration is done by using mb_calibration_from_geodetic_mb to determine the MB model parameters. We will reproduce the option from the preprocessed directory first:

workflow.execute_entity_task(mb_calibration_from_geodetic_mb,
gdir_hef,
informed_threestep=True,   # "informed_threestep" is what OGGM used for 1.6
overwrite_gdir=True)

2023-03-10 19:21:38: oggm.workflow: Execute entity tasks [mb_calibration_from_geodetic_mb] on 1 glaciers

[{'rgi_id': 'RGI60-11.00897',
'bias': 0,
'melt_f': 5.0,
'prcp_fac': 3.5324792824205744,
'temp_bias': 1.758147556769245,
'reference_mb': -1100.3,
'reference_mb_err': 171.8,
'reference_period': '2000-01-01_2020-01-01',
'temp_all_solid': 0.0,
'temp_all_liq': 2.0,
'temp_melt': -1.0},
'baseline_climate_source': 'GSWP3_W5E5'}]

# if you want to know more about that function
# check out the docstring:
# mb_calibration_from_geodetic_mb?


Now how did OGGM come to these parameters?

### First step: data informed precipitation factor#

prcp_fac is chosen from a relation to the average winter precipitation. Fig.1a (below) shows the used relationship between winter precipitation and prcp_fac.

It was calibrated by adapting melt_f and prcp_fac to match the average geodetic and winter MB on around 100 glaciers with both informations available. The found relationship of decreasing prcp_fac for increasing winter precipitation makes sense, as glaciers with already a large winter precipitation should not be corrected with a large multiplicative prcp_fac.

### Second step: data informed temperature bias#

Using this data-informed precipitation factore, with then do a global calibration where the temperature bias (temp_bias) is calibrated, while the melt factor (melt_f) is fixed at 5 kg m-2 day-1 K-1 (default value based on Schuster et al (2023, in review)).

The idea is that if many glaciers within the same grid point need a temperature bias to reach the obeserved MB, this indicates that a systematic correction is necessary (at least for this MB model in particular). In fact, we can plot the median bias required to match MB pbservations using this technique, which gives us the following plot:

The fact that the temp_bias parameter is spatially correlated (many regions are all blue or red) indicate that something in the data needs to be corrected for our model. It is this imformation that we use to inform the next step.

### Third step: local calibration using the data informed precipitation factor and temperature bias as first guesses#

Finally, we now run mb_calibration_from_scalar_mb again for each glacier, as follows:

• use the first guess: melt_f = 5, prcp_fac = data-informed from step 1, temp_bias = data-informed from step 2

• if this doesn’t match (this would be highly unlikely), allow prcp_fac to vary again between 0.8 and 1.2 times the roiginal guess ($$\pm$$20%). This is justified by the fact that the first guess for precipitation is also highly uncertain. If that worked, the calibration stops (XX% of all glaciers worldwide are calibrated this way).

• if the above did not work, allow melt_f to vary again. If that worked, the calibration stops (XX% of all glaciers worldwide are calibrated this way).

• finally, if the above did not work, allow temp_bias to vary again.

This informed step-wise calibration is the option that is used operationally in OGGM>=v1.6 in the preprocessed directories (L>=3) for all glaciers world-wide!

### A look into the global parameter distributions#

w_prcp_array = np.arange(0.5,20,1)
# we basically do here the same as in massbalance.decide_winter_precip_factor(gdir)
a, b = cfg.PARAMS['winter_prcp_fac_ab']
r0, r1 = cfg.PARAMS['prcp_fac_min'], cfg.PARAMS['prcp_fac_max']
prcp_fac = a * np.log(w_prcp_array) + b
# don't allow extremely low/high prcp. factors!!!
prcp_fac_array = utils.clip_array(prcp_fac, r0, r1)
plt.figure(figsize=(12, 5))
plt.subplot(121)
plt.plot(w_prcp_array, prcp_fac_array)
plt.xlabel(r'winter daily mean precipitation' +'\n'+r'(kg m$^{-2}$ day$^{-1}$)')
plt.ylabel('precipitation factor (prcp_fac)'); plt.title('Fig. 1a');
plt.subplot(122)
temp_b_precalib = utils.get_temp_bias_dataframe()['median_temp_bias_w_err_grouped']
plt.hist(temp_b_precalib);
plt.xlabel('temp_bias (°C)'); plt.ylabel('gridpoint frequency'); plt.title('Fig. 1b');


TODO: here we could also add the global distrib of melt_f and prcp_f btw

## Take home points#

• We illustrated how a mass-balance (MB) model of a glacier can be calibrated in OGGM.

• We can use different observational data for calibration:

• calibrating to geodetic observations using different time periods, to in-situ direct glaciological observations from the WGMS (if available) or to other custom MB data.

• There exist different ways of calibrating to the average observational data:

• You can use the same option (informed_threestep) as done operationally for the preprocessed directories (L>=3) on all glaciers world-wide

• you can calibrate the melt_f, and have the prcp_fac and temp_bias fixed. If the calibration does not work, the temp_bias can be varied aswell.

• you can also calibrate instead on prcp_fac or temp_bias and fix the other parameters.

• However, we showed that the parameter combination choice has an influence on other estimates than the average MB (which then influences also future projections).

• As user of OGGM, you might just use the calibrated MB model parameters from the preprocessed directories. Nevertheless, it is good to be aware of the overparameterisation problem. In addition, if you want to include uncertainties of the MB model calibration, you could include additional experiments that use another calibration option and create your own preprocessed directories with these options. With more available observational data or improved climate data, you might also be able to use better ways to calibrate the MB model parameters.