Dynamic spinup for past simulations

New in version 1.5.3!

One of the issues facing glacier modelers is that the past state of glaciers is unknown.

In fact, reconstructing the past states of glaciers turns out to be a really difficult task, and OGGM team member Julia Eis spent her entire PhD thesis doing just that (see this paper for the theory and this one for the application).

However, running simulations in the recent past can be quite useful for model validation. Another very useful use case is to release a strong assumption in OGGM default settings: that glaciers are in dynamical equilibrium at the glacier outline date (an assumption required for the ice thickness inversion).

In recent PRs (GH1342, GH1232 and GH1361), we have released a new run task in OGGM that helps with these issues.

High level explanation of how the method works:

The dynamic spinup iteratively tries to find an initial glacier state in the past (t_start) that match the area or the volume at the glacier outline date (t_end) under the historical climate forcing. To create changing initial glacier states, a constant mass balance is used (mb_spinup) together with a changing temperature bias. This mass balance defines an average mass balance between t_start and t_end. Afterwards, for the actual run, the method uses the OGGM default mass balance model (mb_historical). One iteration of the dynamic spinup looks like this:

  • Apply a temperature bias to mb_spinup.

  • Start with the glacier geometry after the ice thickness inversion and let the model evolve using mb_spinup for the same amount of time as the following historical run (= t_end - t_start).

  • The resulting glacier geometry is defined as the initial glacier state at t_start.

  • Now, let the glacier evolve using mb_historical from t_start until t_end

  • Get the model area or volume of the glacier at t_end.

  • Compare the model value to the reference value we want to meet.

  • If the difference is inside a given precision (precision_percent) stop the procedure and save the glacier evolution of this run (between t_start and t_end).

  • If the difference is outside a given precision, change the temperature bias for mb_spinup and start over again.

To start, the first guess temperature bias is used (first_guess_t_bias) and evaluated. If by chance, the mismatch is fitting, the algorithm stops already. Otherwise, the second guess depends on the calculated mismatch. If the first resulting area or volume is smaller (larger) than the searched one, the second temperature bias will be colder (warmer). Because a colder (warmer) temperature leads to a larger (smaller) initial glacier state at t_start. If still unsuccessful, the previous value pairs (temperature bias, mismatch in per cent) are used to determine the next guess. For this, a stepwise linear function is fitted to these pairs and afterwards, the mismatch is set to 0 to get the following temperature bias (this method is similar to the one described in Zekollari et al. 2019 Appendix A). Moreover, a maximum step length between two guesses is defined with t_bias_max_step_length as too large step sizes could lead to one of the two main problems explained below.

Two main problems why the dynamic spinup could not work:

The resulting glacier is too large, even when we start from a completely ice-free initial glacier state. Or the resulting glacier is too small and the algorithm tries to let the glacier grow outside the domain boundaries. The second problem could be met with a larger domain boundary (if available and enough memory resources available). But for the first problem, no such solution is possible. To still try to cope with these two issues the dynamic spinup function tries to shorten the spinup period two times. The idea is with a shorter evolution time the glacier has less amount of time to grow/shrink and maybe do not encounter one of the two extremes. The minimum length of the spinup can be defined with min_spinup_period (default 10 years). If nothing helps and there is still an error (and if ignore_errors = True) the glacier geometry without a dynamic spinup is saved (= the glacier geometry after the ice thickness inversion).

¡Careful! This method is still a bit experimental. A few things to remember when using it:

  • because of how the method works, glacier area and volume are not strictly the same at the glacier inventory date than after the inversion. We have checked that the errors are random: if you simulate many glaciers this should be fine (global plots on this point will follow soon).

  • because of the non-unique solutions to this optimization problem, the estimated past states become increasingly uncertain the further back in time we go. We recommend a spinup of about 10-30 years.

import matplotlib.pyplot as plt
import xarray as xr
import numpy as np
import pandas as pd
import seaborn as sns
# Make pretty plots
from oggm import cfg, utils, workflow, tasks, graphics
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGMSpinup')
2022-04-03 16:02:28: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2022-04-03 16:02:28: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2022-04-03 16:02:28: oggm.cfg: Multiprocessing: using all available processors (N=2)

Define the glacier we will play with

# Hintereisferner
rgi_id = 'RGI60-11.00897'

# Baltoro
# rgi_id = 'RGI60-14.06794'

Preparing the glacier data

This can take up to a few minutes on the first call because of the download of the required data:

# We use a recent gdir setting, calibated on a glacier per glacier basis
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/ERA5/elev_bands/qc0/pcp1.6/match_geod_pergla/'

# These parameters correspond to the settings of the base url
cfg.PARAMS['prcp_scaling_factor'] = 1.6
cfg.PARAMS['climate_qc_months'] = 0
cfg.PARAMS['hydro_month_nh'] = 1
cfg.PARAMS['hydro_month_sh'] = 1
2022-04-03 16:02:29: oggm.cfg: PARAMS['prcp_scaling_factor'] changed from `2.5` to `1.6`.
2022-04-03 16:02:29: oggm.cfg: PARAMS['climate_qc_months'] changed from `3` to `0`.
2022-04-03 16:02:29: oggm.cfg: PARAMS['hydro_month_nh'] changed from `10` to `1`.
2022-04-03 16:02:29: oggm.cfg: PARAMS['hydro_month_sh'] changed from `4` to `1`.
# We use a relatively large border value to allow the glacier to grow during spinup
gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=5, prepro_border=80, prepro_base_url=base_url)[0]
2022-04-03 16:02:29: oggm.workflow: init_glacier_directories from prepro level 5 on 1 glaciers.
2022-04-03 16:02:29: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

A simple spinup run and projection

# We will "reconstruct" a possible glacier evolution from this year onwards
spinup_start_yr = 1979

The simulation will now consist of two runs:

  • the spinup, from the spinup year to the end of the RGI year (2003)

  • the recent past, from the RGI date to present

We have to use some of the OGGM semantics to get this done:

# New! The dynamic spinup
                         spinup_start_yr=spinup_start_yr,  # When to start the spinup
                         minimise_for='area',  # what target to match at the RGI date
                         output_filesuffix='_spinup_dynamic',  # Where to write the output - this is needed to stitch the runs together afterwards

# This is the same as before - the only difference is that we use the end of the spinup as starting geometry
                            init_model_filesuffix='_spinup_dynamic',  # Which initial geometry to use (from the spinup here: default is from the inversion)
                            output_filesuffix='_hist_spinup_dynamic',  # Where to write the output

# New: stich the output together for analysis
ds_spinup_dynamic = tasks.merge_consecutive_run_outputs(gdir,
                                                        input_filesuffix_1='_spinup_dynamic',  # File 1
                                                        input_filesuffix_2='_hist_spinup_dynamic',  # File 2
                                                        output_filesuffix='_merged_spinup_dynamic',  # Output file (optional)

For comparison, let’s also do a run starting from the inversion geometry. We also use the “fixed geometry spinup” option which is a simpler way to deal with the issue:

# The previous (simple) way of doing this: fixed geometry
                            fixed_geometry_spinup_yr=spinup_start_yr,  # Start the run at the RGI date but retro-actively correct the data with fixed geometry
                            output_filesuffix='_hist_fixed_geom',  # where to write the output
# Read the output
with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_hist_fixed_geom')) as ds:
    ds_hist = ds.load()
y0 = gdir.rgi_date + 1

f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 5))
ds_spinup_dynamic.volume_m3.plot(ax=ax1, label='Dynamical spinup');
ds_hist.volume_m3.plot(ax=ax1, label='Fixed geometry spinup');
ax1.set_title('Volume'); ax1.legend();
ax1.scatter(y0, ds_hist.sel(time=y0).volume_m3, c='C3')
ax2.scatter(y0, ds_hist.sel(time=y0).area_m2, c='C3')
ax3.scatter(y0, ds_hist.sel(time=y0).length_m, c='C3')

In this example the two simulations are very close. But this is actually by chance - other settings will lead to other results, for example starting at another date or optimizing for volume instead of area match. The orange curve shows how the fixed geometry spinup ignores area changes before the RGI date. The short advance after the RGI date is the “initialisation shock”, the issue that the dynamical spinup wants to address.

TODO: add more examples here, including one which doesn’t work that well.

Dynamical spinup and hydrological output

Fortunately, it is also possible to run the spinup model with additional hydrological output:

cfg.PARAMS['store_model_geometry'] = True  # This is required for the hydro runs

# Dynamic spinup
                     run_task=tasks.run_dynamic_spinup,  # The task to run with hydro
                     spinup_start_yr=spinup_start_yr,  # When to start the spinup
                     minimise_for='area',  # what target to match at the RGI date
                     store_monthly_hydro=True,  # compute monthly hydro diagnostics
                     ref_area_from_y0=True,  # Even if the glacier may grow, keep the reference area as the year 0 of the simulation
                     output_filesuffix='_spinup_dynamic_hydro',  # Where to write the output - this is needed to stitch the runs together afterwards

# Recent past 
                     run_task=tasks.run_from_climate_data,  # The task to run with hydro
                     init_model_filesuffix='_spinup_dynamic_hydro',  # Use the spinup as initial conditions!
                     store_monthly_hydro=True,  # compute monthly hydro diagnostics
                     ref_geometry_filesuffix='_spinup_dynamic_hydro',  # Also use the spinup as reference hudrological area!
                     ref_area_from_y0=True,  # Even if the glacier may grow, keep the reference area as the year 0 of the simulation
                     output_filesuffix='_recent_hydro',  # Where to write the output - this is needed to stitch the runs together afterwards

# New: stich the output together for analysis
ds = tasks.merge_consecutive_run_outputs(gdir,
                                         input_filesuffix_1='_spinup_dynamic_hydro',  # File 1
                                         input_filesuffix_2='_recent_hydro',  # File 2

ds = ds.isel(time=slice(0, -1))  # The last timestep is incomplete for hydro (not started)
2022-04-03 16:02:31: oggm.cfg: PARAMS['store_model_geometry'] changed from `False` to `True`.

Annual runoff

The new hydrological variables are now available. Let’s make a pandas DataFrame of all “1D” (annual) variables:

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

The hydrological variables are computed on the largest possible area that was covered by glacier ice in the simulation. This is equivalent to the runoff that would be measured at a fixed-gauge station at the glacier terminus. The total annual runoff is:

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

It consists of the following components:

  • melt off-glacier: snow melt on areas that are now glacier free (i.e. 0 in the year of largest glacier extent, in this example at the start of the simulation)

  • melt on-glacier: ice + seasonal snow melt on the glacier

  • liquid precipitaton on- and off-glacier (the latter being zero at the year of largest glacial extent, in this example at start of the simulation)

f, ax = plt.subplots(figsize=(10, 6));
df_runoff.plot.area(ax=ax, color=sns.color_palette("rocket")); plt.xlabel('Years'); plt.ylabel('Runoff (Mt)'); plt.title(rgi_id);

As the glacier retreats, total runoff decreases as a result of the decreasing glacier contribution.

Monthly runoff

The “2D” variables contain the same data but at monthly resolution, with the dimension (time, month). For example, runoff can be computed the same way:

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

See hydrological_output for more info on how to interpret the hydrological output!

What’s next?