Temperature index models#

Goals of this notebook:

  • Gain a basic understanding of temperature index models

  • Implement OGGM’s temperature index model for a glacier of interest

This version of the notebook works for OGGM versions before 1.6. We will keep this notebook here for a while longer (e.g.: for classroom.oggm.org), and we will replace it with an updated notebook at a later stage. If you are running OGGM 1.6, you should ignore this notebook. 
# Plotting libraries and plot style
import matplotlib.pyplot as plt

import seaborn as sns
sns.set_context('notebook')
sns.set_style('ticks')

import numpy as np
import oggm
from oggm import utils, cfg, workflow, graphics

cfg.initialize()

2023-08-27 22:06:11: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-08-27 22:06:11: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-08-27 22:06:11: oggm.cfg: Multiprocessing: using all available processors (N=2)

Some settings:

# define a temporary directory to store the required data to
cfg.PATHS['working_dir'] = utils.gettempdir('ti_model')

# set the size of the local glacier map: number of grid points outside the
# glacier boundaries
# increasing this parameter will (significantly!) increase the amount of data
# that needs to be downloaded
cfg.PARAMS['border'] = 10

2023-08-27 22:06:13: oggm.cfg: PARAMS['border'] changed from `80` to `10`.

Background#

Glacier melt significantly influences catchment hydrology. Hence, it is useful to have accurate predictions of runoff from glacierized areas. Generally, there are two classes of melt models:

  • energy balance models

  • temperature index models

Energy balance models are physical models quantifying melt as the residual of the energy balance equation. These models require measurements of net radiation, wind speed, temperature and surface properties to predict melt. On a glacier, spatially well resolved measurements are demanding and hard to maintain. Hence, a simpler model, the temperature index model, is the most common approach to model glacier melt.

Temperature index models assume an empirical relationship between air temperatures and melt rates and are a simplification of the energy balance models. The reasoning is that melt is predominantly influenced by the longwave atmospheric radiation and the sensible heat flux - energy balance components that are highly influenced by air temperature (Hock, 2003). The main reason(s) why temperature index models are commonly used are the wide availability of air temperature measurements and computational efficiency.

Model setup#

The simplest temperature index model relates the amount of ice or snow melt \(M\) (mm) to the sum of positive air temperatures \(T^+\) (\(^\circ\)C) by a proportionality factor \(DDF\), the degree-day factor, for each \(n\) time intervals \(\Delta t\):

\[\sum_i^{n} M = DDF \sum_i^{n} T^+ \Delta t\]

Commonly, \(\Delta t = 1\) day is used - hence the name degree-day factor. However, any other time interval \(\Delta t\), e.g. hourly or monthly, can be used to determine \(DDF\). In practice, the model requires measurements of air temperature and glacier mass balance to estimate \(DDF\) - once calculated, \(DDF\) can be used to predict melt by only measuring air temperature (Hock, 2003). However, this temperature index model, also called degree-day model, is not able to predict glacier surface mass balance.

To model glacier surface mass balance, a more sophisticated temperature index model was developed by Marzeion et al., (2012). The monthly mass balance \(B_i\) at elevation \(z\) is computed as

\[B_i(z) = P_i^{solid}(z) - \mu^* \text{max}(T_i(z) - T_{melt}, 0) - \epsilon\]

where \(P_i^{Solid}\) is the monthly solid precipitation, \(T_i\) the monthly average temperature, \(T_{Melt}\) the monthly average temperature above which ice melt is assumed and \(\epsilon\) the residual. \(\epsilon\) is assumed to be a random error taking account for uncertainties associated with unresolved physical processes. \(\mu^*\) is the temperature sensitivity of the glacier and it depends on many parameters, mostly glacier specific (e.g., avalanches, topographical shading, cloudiness, …).

Degrees of freedom#

Among others, the temperature sensitivity \(\mu^*\), the threshold for melt \(T_{Melt}\) and the implicit threshold for solid precipitation \(T_{Solid}\) are important degrees of freedom of the model - \(T_{Solid}\) is the monthly average temperature below which precipitation is assumed to be solid.

Generally, \(T_{Melt}\) and \(T_{Solid}\) can vary both spatially and temporally on a specific glacier. However, commonly the two thresholds \(T_{Melt}\) and \(T_{Solid}\) are assumed to be constant. \(T_{Melt}\) and \(T_{Solid}\) significantly influence the predicted mass balance \(B\) by determining the months which are taken into account in the calculation.

Both \(T_{Melt}\) and \(T_{Solid}\) can be determined by a physical reasoning: we know that both snow melts and precipitation becomes solid at around 0\(^{\circ}\)C. Hence, the two thresholds \(T_{Melt}\) and \(T_{Solid}\) are within a natural range that depends on the climatological conditions at a specific glacier site.

In OGGM, \(T_{Melt}\) and \(T_{Solid}\) are constants and you can access the default values via the cfg module:

# the default temperature below which solid precipitation is assumed
print('T_solid = {}°C'.format(cfg.PARAMS['temp_all_solid']))
# the default temperature above which melt is assumed to occur
print('T_melt = {}°C'.format(cfg.PARAMS['temp_melt']))

T_solid = 0.0°C
T_melt = -1.0°C

Similarly, you can use your own \(T_{Melt}\) and \(T_{Solid}\) if you feel like it:

# don't run this ...
# cfg.PARAMS['temp_all_solid'] = 100
# cfg.PARAMS['temp_melt'] = - 273.15

The temperature sensitivity \(\mu^*\) is glacier specific and mostly determined using statistical error minimization techniques, e.g. ordinary least squares (OLS). Such statistical techniques are very sensitive to the sample size - a general issue in glaciology is that the sample size of annual mass balance records is poor for many glaciers.

However, assume that a \(100\) year long mass balance record together with temperature and precipitation measurements is available for a specific glacier (this is a best case example and only very few glaciers actually have such long records). OLS will find a statistically significant \(\mu^*\) which you can happily use to model mass balance. But what happens if you only use \(30\) years out of the \(100\) year record? OLS will find another statistically significant \(\mu^*\) that is different from the one determined by the \(100\) year record - and another statistically significant \(\mu^*\) can be found for each reasonable subset of the original \(100\) year record. This implies that \(\mu^*\) is generally a time dependent temperature sensitivity \(\mu^*(t)\).

For this reason, OGGM implements a calibration procedure, introduced by Marzeion et al., (2012), to determine a constant glacier specific \(\mu^*\) out of the time dependent \(\mu^*(t)\) candidates. This calibration is beyond the scope of this notebook and you can read about it in detail here and check out an example implementation in OGGM here.

Implementation in OGGM#

First, we need to define a glacier directory:

# this may take a while
gdir = workflow.init_glacier_directories([utils.demo_glacier_id('hef')], from_prepro_level=3)[0]

---------------------------------------------------------------------------
InvalidParamsError                        Traceback (most recent call last)
Cell In[6], line 2
      1 # this may take a while
----> 2 gdir = workflow.init_glacier_directories([utils.demo_glacier_id('hef')], from_prepro_level=3)[0]

File /usr/local/pyenv/versions/3.10.12/lib/python3.10/site-packages/oggm/workflow.py:373, in init_glacier_directories(rgidf, reset, force, from_prepro_level, prepro_border, prepro_rgi_version, prepro_base_url, from_tar, delete_tar)
    370     reset = utils.query_yes_no('Delete all glacier directories?')
    372 if from_prepro_level:
--> 373     url = utils.get_prepro_base_url(base_url=prepro_base_url,
    374                                     border=prepro_border,
    375                                     prepro_level=from_prepro_level,
    376                                     rgi_version=prepro_rgi_version)
    377     if cfg.PARAMS['has_internet'] and not utils.url_exists(url):
    378         raise InvalidParamsError("base url seems unreachable with these "
    379                                  "parameters: {}".format(url))

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

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

If you want to look at your model domain, you can plot it using:

# graphics.plot_domain(gdir)

In OGGM, the calibrated temperature index model for each glacier is accessible via the PastMassBalance class of the massbalance module:

from oggm.core import massbalance

# this class is the calibrated temperature index model
mb_cal = massbalance.PastMassBalance(gdir)

In this case,

print('the glacier selected is {},'.format(gdir.name))

and its calibrated temperature sensitivity \(\mu^*\) is

print('mu_star = {:2f} mm K^-1 yr^-1.'.format(mb_cal.mu_star))

Similarly, the residual \(\epsilon\) is

print('epsilon = {:2f} mm.'.format(mb_cal.bias))

Climate data#

Per default, the temperature index model is driven by the \(0.5^{\circ} \times 0.5^{\circ}\) gridded global CRU TS climate dataset. These climate data are then downscaled to a higher resolution grid to allow for an elevation-dependent dataset. The climate data at the reference height used to drive the temperature index model and to determine the calibrated \(\mu^*\) of the selected glacier can be accessed via the glacier directory:

fpath = gdir.get_filepath('climate_historical')

print(fpath)

This is the temporary path where OGGM stored its climate data on your machine. You can read the climate data using xarray:

import xarray as xr

climate = xr.open_dataset(fpath)
climate

The climate dataset has two variables, the monthly total precipitation prcp and the monthly average temperature temp. Let’s calculate the mean annual cycle of average temperature and total precipitation,

annual_cycle = climate.groupby('time.month').mean(dim='time')

and plot it, to get an intuitive view on the climate conditions at the selected glacier site.

import calendar

fig, ax = plt.subplots(1, 2, figsize=(16, 9))
ax[0].plot(annual_cycle.month, annual_cycle.temp); ax[1].plot(annual_cycle.month, annual_cycle.prcp);
ax[0].set_title('Average temperature / (°C)'); ax[1].set_title('Total precipitation / (mm)');
for a in ax:
    a.set_xticks(annual_cycle.month.values)
    a.set_xticklabels([calendar.month_abbr[m] for m in annual_cycle.month.values])

Reference mass balance data#

OGGM uses in-situ mass balance data from the World Glacier Monitoring Service Fluctuations of Glaciers Database (WGMS FoGD). The Fluctuations of Glaciers (FoG) database contains annual mass-balance records for several hundreds of glaciers worldwide. Currently, 254 mass balance time series are available.

These data are shipped automatically with OGGM and can be accessed via the glacier directory:

# Get the reference mass-balance from the WGMS FoG Database
ref_mb = gdir.get_ref_mb_data()

ref_mb[['ANNUAL_BALANCE']].plot(title='Annual mass balance: {}'.format(gdir.name), legend=False);

Predict mass balance!#

Now, we are set to calculate glacier mass balance using the temperature index model - we have the model parameters \(\mu^*\) and \(\epsilon\), the thresholds for melt and solid precipitation \(T_{Melt}\) and \(T_{Solid}\) and the climate dataset. The last thing we need to define are the heights at which we want to calculate the mass balance. Here, we use glacier flowlines along which the mass balance is computed:

fls = gdir.read_pickle('inversion_flowlines')

We will calculate the specific mass balance in mm w.e. yr\(^{-1}\) for the years where in-situ mass balance data is available:

print(ref_mb.index.values)

The specific mass balance along the given flowlines is computed by

ref_mb['OGGM (calib)'] = mb_cal.get_specific_mb(fls=fls, year=ref_mb.index.values)

For this calculation we assumed an average ice density of

print('rho_ice = {} kg m^-3.'.format(cfg.PARAMS['ice_density']))

Now, we can compare the actual in-situ mass balance with the modelled mass balance:

fig, ax = plt.subplots(1, 1)
ax.plot(ref_mb['ANNUAL_BALANCE'], label='Observed')
ax.plot(ref_mb['OGGM (calib)'], label='Modelled')
ax.set_ylabel('Specific mass balance / (mm w.e. y$^{-1}$)')
ax.legend(frameon=False);

Does not look too bad, does it? To assess model performance, it is helpful to plot the data in a scatter plot:

fig, ax = plt.subplots(1, 1)
ax.plot(ref_mb['ANNUAL_BALANCE'], ref_mb['OGGM (calib)'], 'ok');
ax.plot(ref_mb['ANNUAL_BALANCE'], ref_mb['ANNUAL_BALANCE'], '-r');
ax.set_xlim(-3000, 2000)
ax.set_ylim(-3000, 2000)
ax.set_xlabel('Observed');
ax.set_ylabel('OGGM (calib)');

If the points were aligned along the red line, the model would perfectly predict mass balance. Generally, the model overestimates mass balance in magnitude - the scatter plot shows a steeper slope than the 1 to 1 red line. This is due to specific mass balance beeing dependent not only on the climate but also on the glacier surface area.

OGGM computes the specific mass balance as a glacier area-weighted average using a constant glacier geometry fixed at the Randolph Glacier Inventory date, e.g. \(2003\) for most glaciers in the European Alps. Glacier geometry is itself a function of climate and may change significantly over time. Hence, assuming a constant glacier geometry over a time period of different climatic conditions can result in a systematic model bias:

bias = ref_mb['OGGM (calib)'] - ref_mb['ANNUAL_BALANCE']

fig, ax = plt.subplots(1, 1)
ax.plot(bias);

The bias is positive at the beginning of the in-situ measurements and shows a negative trend. When keeping the glacier area constant, a positive (negative) bias means, that the calibrated temperature sensitivity \(\mu^*\) of the glacier is too low (high) during time periods of colder (warmer) climates. You can find a simple experiment about the sensitivity of the specific mass balance on climate change and glacier surface area in this blog post.

Take home points#

  • There are two different types of melt models: the energy balance model and the temperature index model

  • The temperature index model is the computationally efficient simplification of the energy balance model

  • Temperature index models assume an empirical relationship between air temperature and melt rates

  • Temperature index models can be extended to model glacier mass balance by adding solid precipitation as a model parameter

  • The model outcome is significantly influenced by the choice of \(T_{Melt}\) and \(T_{Solid}\)

  • The temperature sensitivity of a glacier is not constant in time \(\mu^* = \mu^*(t)\)

  • The specific mass balance is a function of the climate and the glacier surface area

References#

What’s next?#

Back to the table of contents