Plotting the OGGM surface mass-balance, the ELA and AAR#

A few recipes to help you to plot and analyse the mass-balance computed by OGGM.

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

import oggm
from oggm import cfg, utils, workflow, tasks, graphics, global_tasks
from oggm.core import massbalance, flowline

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

2023-03-07 12:45:40: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-03-07 12:45:40: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-07 12:45:40: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-03-07 12:45:40: oggm.cfg: PARAMS['border'] changed from `40` to `80`.
2023-03-07 12:45:40: oggm.cfg: PARAMS['store_model_geometry'] changed from `False` to `True`.

We use Hintereisferner:

# (The second glacier here we will only use in the second part of this notebook.)
gdirs = workflow.init_glacier_directories(['RGI60-11.00897', 'RGI60-11.00001'], from_prepro_level=3)
gdir = gdirs[0]
tasks.init_present_time_glacier(gdir)

2023-03-07 12:45:41: oggm.workflow: init_glacier_directories from prepro level 3 on 2 glaciers.
2023-03-07 12:45:41: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 2 glaciers

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[3], line 2
      1 # (The second glacier here we will only use in the second part of this notebook.)
----> 2 gdirs = workflow.init_glacier_directories(['RGI60-11.00897', 'RGI60-11.00001'], from_prepro_level=3)
      3 gdir = gdirs[0]
      4 tasks.init_present_time_glacier(gdir)

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

With a static geometry#

That’s by far the easiest - that is how we mostly test and develop mass-balance models ourselves.

The very easy way#

# Get the calibrated mass-balance model - the default is to use OGGM's "PastMassBalanceModel"
mbmod = massbalance.MultipleFlowlineMassBalance(gdir)

# Compute the specific MB for this glacier
fls = gdir.read_pickle('model_flowlines')
years = np.arange(1950, 2019)
mb_ts = mbmod.get_specific_mb(fls=fls, year=years)

plt.plot(years, mb_ts);
plt.ylabel('Specific MB (mm w.e.)');
plt.xlabel('Year');

In this particular case (Hintereisferner), we can also plot the measured data:

# Get the reference mass-balance from the WGMS
ref_df = gdir.get_ref_mb_data()
# Get OGGM MB at the same dates
ref_df['OGGM'] = mbmod.get_specific_mb(fls=fls, year=ref_df.index.values)
# Plot
plt.plot(ref_df.index, ref_df.OGGM, label='OGGM');
plt.plot(ref_df.index, ref_df.ANNUAL_BALANCE, label='WGMS');
plt.ylabel('Specific MB (mm w.e.)'); plt.legend();
plt.xlabel('Year');

Mass-balance profiles on flowlines#

Let’s start by noting that I used the MultipleFlowlineMassBalance model here. This is what OGGM uses for its runs, because we allow different model parameters for different flowlines. In many cases we don’t need this, but sometimes we do. Let’s see if we have different temperature sensititivites for Hintereisferner or not:

for flmod in mbmod.flowline_mb_models:
    print(flmod.mu_star)

So here we have the same model parameters for each flowline (that’s frequently the case). To get the mass-balance as a function of height we have several possibilities:

# Where the flowlines are:
heights, widths, mb = mbmod.get_annual_mb_on_flowlines(fls, year=2000)
# units kg ice per second to mm w.e. per year:
mb = mb * cfg.SEC_IN_YEAR * cfg.PARAMS['ice_density'] 

# Plot
plt.plot(mb, heights, '*', label='2000');

# Another year:
heights, widths, mb = mbmod.get_annual_mb_on_flowlines(fls, year=2001)
# units kg ice per second to mm w.e. per year:
mb = mb * cfg.SEC_IN_YEAR * cfg.PARAMS['ice_density'] 
# Plot
plt.plot(mb, heights, '*', label='2001');
plt.ylabel('Elevation (m a.s.l.)'); plt.xlabel('MB (mm w.e. yr$^{-1}$)'); plt.legend();

Mass-balance at any height#

For this, we can either “pick” one of the single flowline mass-balance models or use the multiple one (it’s the same):

heights = np.arange(2500, 4000)
# Method 1
mbmod1 = mbmod.flowline_mb_models[0]
mb = mbmod1.get_annual_mb(heights, year=2001) * cfg.SEC_IN_YEAR * cfg.PARAMS['ice_density'] 
# Method 2
mb2 = mbmod.get_annual_mb(heights, year=2001, fl_id=0) * cfg.SEC_IN_YEAR * cfg.PARAMS['ice_density'] 
np.testing.assert_allclose(mb, mb2)
# Plot
plt.plot(mb, heights, label='2001');
plt.ylabel('Elevation (m a.s.l.)'); plt.xlabel('MB (mm w.e. yr$^{-1}$)'); plt.legend();

But again, for certain glaciers the mass-balance might change from one flowline to another!

Mass-balance during a transient run + bonus info (timestamps in OGGM)#

Specific (glacier-wide) MB#

At the time of writing, we only store the glacier volume, which is roughly equivalent to the mass-balance (up to numerical accuracy in the dynamical model).

# We want to use these flowline files for the AAR (see later)
cfg.PARAMS['store_fl_diagnostics'] = True
# Run from outline inventory year (2003 for HEF) to 2017 (end of CRU data in hydro year convention, see below)
tasks.run_from_climate_data(gdir, ye=2020);

ds_diag = xr.open_dataset(gdir.get_filepath('model_diagnostics'))
ds_diag.volume_m3.plot(); plt.title('Volume evolution of Hintereisferner');

Important: timestamps in OGGM#

This can be a bit of a headache, so let’s spend some time explaining this. OGGM was developed to follow the traditional “hydrological year” convention. So in the OGGM output files (and in the model internals), we use the same convention (history will tell if it was a good idea, I sometimes regret it).

We use what we call a “floating year” time stamp, which can be converted to real dates easily in the model:

from oggm.utils import floatyear_to_date, hydrodate_to_calendardate
print('Hydro dates')
print(floatyear_to_date(2003.0))
print(floatyear_to_date(2003.99))
print('Calendar dates')
print(hydrodate_to_calendardate(*floatyear_to_date(2003.0), start_month=10))
print(hydrodate_to_calendardate(*floatyear_to_date(2003.99), start_month=10))

The timestamp used in the OGGM output files uses the same convention:

print(ds_diag.time.values[0], ds_diag.calendar_year.values[0], ds_diag.calendar_month.values[0])

So the simulation started in year 2004.0 (as of the OGGM update from 19.05.2020), which corresponds to the real date October 2003. This is in accordance with the RGI inventory date (probably based on satellite imagery taken in the summer 2003):

gdir.rgi_date

Similarly, the simulation ends at the last possible date for which we have climate data for (i.e. the simulation stops at the end of September 2019, which translates to the beginning of hydrological year 2020):

print(ds_diag.time.values[-1], ds_diag.calendar_year.values[-1], ds_diag.calendar_month.values[-1])

Therefore, you have to be careful about how you interpret these data or report them to others.

Specific MB from volume change time series#

We compute the specific MB from the delta volume divided by the glacier area, and compare it to observations and to the mass-balance computed with a fixed geometry:

# Specific MB over time (i.e. with changing geometry feedbacks)
smb = (ds_diag.volume_m3.values[1:] - ds_diag.volume_m3.values[:-1]) / ds_diag.area_m2.values[1:]
smb = smb * cfg.PARAMS['ice_density']  # in mm
plt.plot(ds_diag.time[:-1], smb, label='OGGM (dynamics)'); 
# The SMB from WGMS and fixed geometry we already have
plt.plot(ref_df.loc[2004:].index, ref_df.loc[2004:].OGGM, label='OGGM (fixed geom)');
plt.plot(ref_df.loc[2004:].index, ref_df.loc[2004:].ANNUAL_BALANCE, label='WGMS');
plt.legend();

A few comments:

  • the systematic difference between WGMS and OGGM is due to the calibration of the model, which is done with fixed geometry (the bias is zero over the entire calibration period but changes sign over the 50 year period, see plot above)

  • the difference between OGGM dynamics and fixed geometry (so small that they are barely visible for this short time period) is due to:

    • the changes in geometry during the simulation time (i.e. this difference grows with time)

    • melt at the tongue which might be larger than the glacier thickness is not accounted for in the fixed geometry data (these are assumed to be small)

    • numerical differences in the dynamical model (these are assumed to be small as well)

MB profiles#

At time of writing, we do not store the mass-balance profiles during a transient run, although we could (and should). To get the profiles with a dynamically changing geometry, there is no other way than either get the mass-balance during the simulation, or to fetch the model output and re-compute the mass-balance retroactively.

Note: the MB profiles themselves of course are not really affected by the changing geometry - what changes from one year to another is the altitude at which the model will compute the MB. In other terms, if you are interested in the mass-balance profiles you are better of to use the fixed geometries approaches explained at the beginning of the notebook.

Retrospective MB from model output#

# This reads the model output file and makes it look like an OGGM flowline model
dyn_mod = flowline.FileModel(gdir.get_filepath('model_geometry'))
# Get the calibrated mass-balance model
mbmod = massbalance.MultipleFlowlineMassBalance(gdir)
# We can loop over the years and read the flowlines for each year
for year in range(2004, 2017):
    dyn_mod.run_until(year)
    # Compute the SMB from the mass balance model and the flowlines heights at each year
    h = []
    smb = []
    for fl_id, fl in enumerate(dyn_mod.fls):
        h = np.append(h, fl.surface_h)
        mb = mbmod.get_annual_mb(fl.surface_h, year=year, fl_id=fl_id) 
        smb = np.append(smb, mb * cfg.SEC_IN_YEAR * cfg.PARAMS['ice_density'] )
    plt.plot(smb, h, '.', label=year)
plt.legend();

Compute the MB during the simulation#

# TODO - add the same, but based on the new flowline outputs (like for AAR below)

Equilibrium Line Altitude#

The second part of this notebook shows how you can compute the Equilbrium Line Altitude (ELA, the altitude at which the MB is equal to 0) with OGGM.

As the ELA in OGGM only depends on the mass balance model itself (not on glacier geometry), there is no need to do a full model run to collect these values. This is also the reason why the ELA is no longer a part of the model run output, but is instead a diagnostic variable that can be computed separately since OGGM v1.6.

Compute and compile the ELA#

global_tasks.compile_ela(gdirs, ys=2000, ye=2019);

Per default, the file will be stored in the working directory and named ELA.hdf (you can choose to save the data in a csv file instead by setting the csv keyword to True). The location where the file is saved can be changed by specifying the path keyword.

Let’s read the file and plot it:

# Read the file
ela_df = pd.read_hdf(os.path.join(cfg.PATHS['working_dir'], 'ELA.hdf'))

# Plot it
ela_df.plot();
plt.xlabel('year'); plt.ylabel('ELA [m]');

When looking at many glaciers it might be useful to also plot the mean or average of the data. Let’s plot the mean ELA and the glacier area weighted ELA average:

areas = [gd.rgi_area_km2 for gd in gdirs]

avg = ela_df.mean(axis=1).to_frame(name='mean')
avg['weighted average'] = np.average(ela_df, axis=1, weights=areas)

avg.plot();
plt.xlabel('year'); plt.ylabel('ELA [m]');

How would to ELA look like if the climate is 1° warmer? Lets have a look at one of the glaciers:

global_tasks.compile_ela(gdirs, ys=2000, ye=2019, temperature_bias=1, filesuffix='_t1')
ela_df_t1 = pd.read_hdf(os.path.join(cfg.PATHS['working_dir'], 'ELA_t1.hdf'))

plt.plot(ela_df.mean(axis=1), label='default climate')
plt.plot(ela_df_t1.mean(axis=1), label=('default climate +1℃'))
plt.xlabel('year CE'); plt.ylabel('ELA [m]'); plt.legend();

By using the precipitation_factor keyword, you can change the precipitation. Feel free to try that by yourself!

Lets look at a longer timeseries, for one glacier only:

global_tasks.compile_ela(gdirs[0], ys=1902, ye=2019, filesuffix='_1901_2019', csv=True)
ela_df_long = pd.read_csv(os.path.join(cfg.PATHS['working_dir'], 'ELA_1901_2019.csv'), index_col=0)

The ELA can have a high year-to-year variability. Therefore we plot in addition to the regular ELA timeseries, the 5-year moving average.

ax = ela_df_long.plot();
ela_df_long.rolling(5).mean().plot(ax=ax, lw=2, color='k', label='')
plt.xlabel('year CE'); plt.ylabel('ELA [m]'); ax.legend(["Annual", "5-yr average"]);

In case you’re only intrested in specific years there is an option to just compute the ELA in those years. There is actually no need to save this data. Therefore we now just compute the ELA.

yrs = [2000, 1950, 1980, 1960]
ela_yrs = massbalance.compute_ela(gdirs[0], years=yrs)
ela_yrs

See in this case you can use any year as long as you have climate available for that year. However for plotting purposes it might be worth to sort the data, otherwise the following happens ;)

ela_yrs.plot();

ela_yrs.sort_index(ascending=True).plot();
plt.xlabel('year CE'); plt.ylabel('ELA [m]');

By now we have addressed most of the keyword agruments of the compile_ela and compute_ela functions. To use different climate than the default climate, you can make use of the climate_filename and climate_input_filesuffix keywords. In case you’re not familiar with those yet, please check out the run_with_gcm notebook.

Accumulation Area Ratio#

Like the ELA, the Accumulation Area Ratio (AAR) is a diagostic variable. The AAR is currently not an output variable in OGGM, but it can easily be computed. It is even a part of the glacier simulator and its documentation. Below we give an example of how it can be computed after a model run. From the ELA computation, we already know the height where the mass balance is equal to zero. Now we need to compute the area of the glacier that is located above that line. Let’s start with a static geometry glacier with multiple flowlines.

We first search the points along the flowlines that are above the ELA in a given year. Lets use Hintereisferner.

We use a True, False array saying which part is above the ELA and sum the corresponding area at each timestep:

rgi_id = 'RGI60-11.00897'

tot_area = 0
tot_area_above = 0
for fl in fls:
    # The area is constant
    tot_area += fl.area_m2
    # The area above a certain limit is not - we make a (nx, nt) array
    # This is a little bit of dimensional juggling (broadcasting)
    is_above = (fl.surface_h[:, np.newaxis, ] > ela_df[rgi_id].values)
    # Here too - but the time dim is constant
    area = fl.bin_area_m2[:, np.newaxis] * np.ones(len(ela_df))
    tot_area_above += (area * is_above).sum(axis=0)

# Write the output
fixed_aar = pd.DataFrame(index=ela_df.index)
fixed_aar['AAR'] = tot_area_above / tot_area

# Plot ELA and AAR on the same plot
fixed_aar['ELA'] = ela_df[rgi_id]
fixed_aar.plot(secondary_y='ELA');

Unsurprisingly, the AAR is strongly correlated to the ELA, but it is not fully linear due to the area distribution of the glacier:

fixed_aar.plot.scatter(x='AAR', y='ELA');

Now that we have looked at the static case, let’s to the same for glaciers with changing surface height:

# We start by reading the flowline geometry diagnostic files
tot_area = 0
tot_area_above = 0

ela_ds = xr.DataArray(ela_df[rgi_id], dims=('time'))

for fn in np.arange(len(fls)):
    with xr.open_dataset(gdir.get_filepath('fl_diagnostics'), group='fl_' + str(fn)) as ds:
        ds = ds.load()
    tot_area += ds.area_m2.sum(dim='dis_along_flowline')

    # Compute the surface elevation
    surface_h = ds.thickness_m + ds.bed_h

    # In xarray things are much easier - the broadcasting is automatic
    is_above = surface_h > ela_ds
    area = is_above * ds.area_m2
    tot_area_above += (area * is_above).sum(dim='dis_along_flowline')

# Write the output
dyna_aar = pd.DataFrame(index=tot_area_above.time)
dyna_aar['AAR_dyna'] = tot_area_above / tot_area
dyna_aar['AAR_fixed'] = fixed_aar['AAR']

dyna_aar.plot(); plt.ylabel('AAR');

You can see in the plot that the difference between the AAR calculated from the dynamic glacier and the one of the static glacier is quite similar. Over time you can expect that this difference becomes larger. Feel free to play around with that.

Keep in mind that this is just an example for computing the AAR. For a different glacier experiment you might need to adjust the code. e.g. if the glacier flowlines would have different temperature sensitivities, you need to take that in to consideration.

Here we computed the AAR from the “model perspective”, where the glacier is located in bins. You can however also choose to do some interpolation between the grid point above and below the ELA.

What’s next?#