Dynamic spinup and dynamic \(\mu^*\) calibration for past simulations#

New in version 1.5.3 resp. version 1.6.0!

One of the issues glacier modellers facing is that the past state of glaciers is unknown. On the other side, more global datasets of past observations get available (e.g. geodetic mass balance), which can help constrain the recent past evolution.

Reconstructing the past states of glaciers on a century timescale 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. Also, more direct observations are available of glacier states of the recent past for constriction (e.g. area, geodetic mass balance). Further, a dynamical initialisation can release strong assumptions in the OGGM default settings: first that glaciers are in dynamical equilibrium at the glacier outline date (an assumption required for the ice thickness inversion) and second that the mass balance temperature sensitivity parameter (\(\mu^*\)) is calibrated towards a geodetic mass balance ignoring a dynamically changing glacier geometry.

In recent PRs (GH1342, GH1232, GH1361 and GH1425) we have released two new run tasks in OGGM which help with this issues:

The run_dynamic_spinup task, by default, aims to find a glacier state before the RGI-date (~10-30 years back) from which the glacier evolves to match the area given by the RGI-outline. Alternatively, it is also possible to use this task to match an observed volume.
The run_dynamic_mu_star_calibration task iteratively searches for a \(\mu^*\) to match the observed geodetic mass balance taking a dynamically changing glacier geometry into account.

And of course, we want to match both things in the same past model run. Therefore by default in each iteration of run_dynamic_mu_star_calibration the run_dynamic_spinup function is included. A more in-depth explanation of the two tasks is provided in the next two chapters, which are followed by an example and a comparison of the different spinup options.

High-level explanation of `run_dynamic_spinup`:#

The dynamic spinup iteratively tries to find an initial glacier state in the past (t_start) that matches 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, stop the procedure and save the glacier evolution of this run.
If the difference is outside a given precision, change the temperature bias for mb_spinup and start over again (how the next guess is found is descript here).

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 are 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 as 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.

High-level explanation of `run_dynamic_mu_star_calibration` :#

This task iteratively searches for a \(\mu^*\) to match a given geodetic mass balance incorporating a dynamic model run. But changing \(\mu^*\) means we need to rerun all model setup steps which incorporate the mass-balance, to have one consistent model initialisation chain. In particular, we need to conduct the bed inversion again (the mass-balance is used in the flux calculation, see here). Therefore one default iteration of the dynamic \(\mu^*\) calibration looks like this:

define a new \(\mu^*\) in the glacier directory
conduct an inversion which calibrates to the consensus volume, still assuming dynamic equilibrium. (a little tricky: by default, before we start with the dynamic \(\mu^*\) calibration the first inversion with calibration on a regional scale was already carried out -> the individual glaciers do not match the consensus volume exactly, but when adding all glaciers of one region the consensus volume is matched. Therefore during this task, the individual glacier is again matched to the volume of the regional assessment and not on an individual basis.)
conduct a dynamic spinup to match the RGI area and let the model evolve up to the end of the geodetic mass balance period
calculate the modelled geodetic mass balance
calculate the difference between modelled value and observation
if the difference is inside the observation uncertainty stop and save the model run (indicated in diagnostics with used_spinup_option = dynamic mu_star calibration (full success))
if the difference is larger, define a new \(\mu^*\) and start over again (how the next guess is found is described here)

If the iterative search is not successful and ignore_errors = True there are several possible outcomes:

First, it is checked if there were some successful runs which improved the mismatch. If so, the best run is saved and it is indicated in the diagnostics with used_spinup_option = dynamic mu_star calibration (part success)
If only the first guess worked this run is saved and indicated in the diagnostics with used_spinup_option = dynamic spinup only
And if everything failed a fixed geometry spinup is conducted and indicated in the diagnostics with used_spinup_option = fixed geometry spinup

If for whatever reason one wants to conduct the dynamic \(\mu^*\) calibration without an inversion or does not want to include a dynamic spinup, reach out to us as there are already options included for this.

The minimization algorithm:#

To start, a first guess of the control variable (temperature bias or \(\mu^*\)) is used and evaluated. If by chance, the mismatch between model and observation is close enough, the algorithm stops already. Otherwise, the second guess depends on the calculated first guess mismatch. For example, if the first resulting area 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 the second guess is still unsuccessful, for all consecutive guesses the previous value pairs (control variable, mismatch) 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 guess (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 as too large step-sizes could easily lead to failing model runs (e.g. see here). Further, the algorithm was adapted independently inside run_dynamic_spinup and run_dynamic_mu_star_calibration to cope with failing model runs individually. Note that this minimization algorithm only works if the underlying relationship between the control variable and the mismatch is strictly monotone.

If someone is interested in how this algorithm works in more detail, here is a conceptual code snippet:

import numpy as np
from scipy import interpolate

def minimisation_algorithm(
        fct_to_minimise,  # the function which should be minimised
        first_guess,  # the first guess of the control variable
        get_second_guess,  # function which gives the second guess depending on
                           # the first mismatch
        precision_to_match,  # we want to be smaller with the mismatch than this
        maxiter,  # how many search iteration should be performed
        ):

    control_guesses = []
    mismatch = []

    # returns also the new control if some error management was going on inside
    # fct_to_minimise, not necessary
    new_mismatch, new_control = fct_to_minimise(first_guess)
    control_guesses.append(new_control)
    mismatch.append(new_mismatch)

    # check if we can stop already
    if abs(mismatch[-1]) < precision_to_match:
            return control_guesses, mismatch

    # how to get the second guess must also be defined by the user
    second_guess = get_second_guess(mismatch[-1])
    new_mismatch, new_control = fct_to_minimise(second_guess)
    control_guesses.append(new_control)
    mismatch.append(new_mismatch)
    if abs(mismatch[-1]) < precision_to_match:
        return control_guesses, mismatch

    # now we start with the partial linear fitting
    while len(control_guesses) < maxiter:
        # sort the (mismatch, control_pairs) from small to large
        sort_index = np.argsort(np.array(mismatch))
        # actual linear fitting
        tck = interpolate.splrep(np.array(mismatch)[sort_index],
                                 np.array(control_guesses)[sort_index],
                                 k=1)  # degree of the fitting, 1 is linear
        # if there are NaNs in the fitting parameters it is an indication that
        # probably the strictly monotone relation is not holding, but also
        # other problems possible (e.g. to close guesses)
        if np.isnan(tck[1]).any():
            # here is the place where you can do further problem dependent
            # error management or just raise an error
            raise RuntimeError('Somthing went wrong during the linear fitting!')

        # evaluate the linear fit where mismatch = 0 to get the next guess
        next_control_guess = float(interpolate.splev(0, tck))

        new_mismatch, new_control = fct_to_minimise(next_control_guess)
        control_guesses.append(new_control)
        mismatch.append(new_mismatch)
        if abs(mismatch[-1]) < precision_to_match:
            return control_guesses, mismatch

    # If we end here it was not possible to find a satisfying mismatch after maxiter(ations)
    raise RuntimeError('Could not find small enough mismatch after the'
                       'maximum iterations!')

Start with glacier example#

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

# Make pretty plots
sns.set_style('ticks')
sns.set_context('notebook')

from oggm import cfg, utils, workflow, tasks, graphics

cfg.initialize(logging_level='WARNING')
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGMSpinup')

2023-03-07 12:25:49: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-03-07 12:25:49: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-07 12:25:49: oggm.cfg: Multiprocessing: using all available processors (N=2)

Define the glaciers we will play with#

# Khumbu and a 'bad' glacier in terms of the dynamic spinup
rgi_ids = ['RGI60-15.03733', 'RGI60-11.00275']

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
cfg.PARAMS['max_mu_star'] = 600

2023-03-07 12:25:49: oggm.cfg: PARAMS['prcp_scaling_factor'] changed from `2.5` to `1.6`.
2023-03-07 12:25:49: oggm.cfg: PARAMS['climate_qc_months'] changed from `3` to `0`.
2023-03-07 12:25:49: oggm.cfg: PARAMS['hydro_month_nh'] changed from `10` to `1`.
2023-03-07 12:25:49: oggm.cfg: PARAMS['hydro_month_sh'] changed from `4` to `1`.
2023-03-07 12:25:49: oggm.cfg: PARAMS['max_mu_star'] changed from `10000.0` to `600`.

# We use a relatively large border value to allow the glacier to grow during spinup
gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_border=160, prepro_base_url=base_url)

---------------------------------------------------------------------------
InvalidParamsError                        Traceback (most recent call last)
Cell In[8], line 2
      1 # We use a relatively large border value to allow the glacier to grow during spinup
----> 2 gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_border=160, prepro_base_url=base_url)

File /usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/workflow.py:493, in init_glacier_directories(rgidf, reset, force, from_prepro_level, prepro_border, prepro_rgi_version, prepro_base_url, from_tar, delete_tar)
    488     url = utils.get_prepro_base_url(base_url=prepro_base_url,
    489                                     border=prepro_border,
    490                                     prepro_level=from_prepro_level,
    491                                     rgi_version=prepro_rgi_version)
    492     if cfg.PARAMS['has_internet'] and not utils.url_exists(url):
--> 493         raise InvalidParamsError("base url seems unreachable with these "
    494                                  "parameters: {}".format(url))
    496 # if reset delete also the log directory
    497 if reset:

InvalidParamsError: base url seems unreachable with these parameters: https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/ERA5/elev_bands/qc0/pcp1.6/match_geod_pergla/RGI62/b_160/L3/

Start by comparing three different spinups#

# For the following experiments we only use Khumbu
gdir = gdirs[0]

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

# here we save the original mu_star for later
mu_star_original = gdir.read_json('local_mustar')['mu_star_glacierwide']

# and save some reference values we can compare to later
area_reference = gdir.rgi_area_m2  # RGI area
volume_reference = tasks.get_inversion_volume(gdir)  # volume after regional calibration to consensus
ref_period = cfg.PARAMS['geodetic_mb_period']  # get the reference period we want to use, 2000 to 2020
df_ref_dmdtda = utils.get_geodetic_mb_dataframe().loc[gdir.rgi_id]  # get the data from Hugonnet 2021
df_ref_dmdtda = df_ref_dmdtda.loc[df_ref_dmdtda['period'] == ref_period]  # only select the desired period
dmdtda_reference = df_ref_dmdtda['dmdtda'].values[0] * 1000  # get the reference dmdtda and convert into kg m-2 yr-1
dmdtda_reference_error = df_ref_dmdtda['err_dmdtda'].values[0] * 1000  # corresponding uncertainty

Let’s start by comparing the fixed_geometry_spinup, the run_dynamic_spinup and the run_dynamic_mu_star_calibration (which contains a dynamic spinup) options.

# ---- First do the fixed geometry spinup ----
tasks.run_from_climate_data(gdir,
                            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()

# ---- Second the dynamic spinup alone, matching area ----
tasks.run_dynamic_spinup(gdir,
                         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_area',  # Where to write the output
                         ye=2020,  # When the simulation should stop
                         );
# Read the output
with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_spinup_dynamic_area')) as ds:
    ds_dynamic_spinup_area = ds.load()
    
# ---- Third the dynamic spinup alone, matching volume ----
tasks.run_dynamic_spinup(gdir,
                         spinup_start_yr=spinup_start_yr,  # When to start the spinup
                         minimise_for='volume',  # what target to match at the RGI date
                         output_filesuffix='_spinup_dynamic_volume',  # Where to write the output
                         ye=2020,  # When the simulation should stop
                         );
# Read the output
with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_spinup_dynamic_volume')) as ds:
    ds_dynamic_spinup_volume = ds.load()

# ---- Fourth the dynamic mu star calibration (includes a dynamic spinup matching area) ----
tasks.run_dynamic_mu_star_calibration(gdir,
                                      ys=spinup_start_yr,  # When to start the spinup
                                      ye=2020,  # When the simulation should stop
                                      output_filesuffix='_dynamic_mu',  # Where to write the output
                                     );

with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_dynamic_mu')) as ds:
    ds_dynamic_mu = ds.load()

# Now make a plot for comparision
y0 = gdir.rgi_date + 1

f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 5))

ds_hist.volume_m3.plot(ax=ax1, label='Fixed geometry spinup');
ds_dynamic_mu.volume_m3.plot(ax=ax1, label='Dynamical mu* calibration');
ds_dynamic_spinup_area.volume_m3.plot(ax=ax1, label='Dynamical spinup match area');
ds_dynamic_spinup_volume.volume_m3.plot(ax=ax1, label='Dynamical spinup match volume');
ax1.set_title('Volume');
ax1.scatter(y0, volume_reference, c='C3', label='Reference values');
ax1.legend();

ds_hist.area_m2.plot(ax=ax2);
ds_dynamic_mu.area_m2.plot(ax=ax2);
ds_dynamic_spinup_area.area_m2.plot(ax=ax2);
ds_dynamic_spinup_volume.area_m2.plot(ax=ax2);
ax2.set_title('Area');
ax2.scatter(y0, area_reference, c='C3')

ds_hist.length_m.plot(ax=ax3);
ds_dynamic_mu.length_m.plot(ax=ax3);
ds_dynamic_spinup_area.length_m.plot(ax=ax3);
ds_dynamic_spinup_volume.length_m.plot(ax=ax3);
ax3.set_title('Length');
ax3.scatter(y0, ds_hist.sel(time=y0).length_m, c='C3')

plt.tight_layout()
plt.show();

# and print out the modeled geodetic mass balances for comparision
def get_dmdtda(ds):
    yr0_ref_mb, yr1_ref_mb = ref_period.split('_')
    yr0_ref_mb = int(yr0_ref_mb.split('-')[0])
    yr1_ref_mb = int(yr1_ref_mb.split('-')[0])

    return ((ds.volume_m3.loc[yr1_ref_mb].values -
             ds.volume_m3.loc[yr0_ref_mb].values) /
            gdir.rgi_area_m2 /
            (yr1_ref_mb - yr0_ref_mb) *
            cfg.PARAMS['ice_density'])

print(f'Reference dmdtda 2000 to 2020 (Hugonnet 2021):     {dmdtda_reference:.2f} +/- {dmdtda_reference_error:6.2f} kg m-2 yr-1')
print(f'Fixed geometry spinup dmdtda 2000 to 2020:         {get_dmdtda(ds_hist):.2f}            kg m-2 yr-1')
print(f'Dynamical spinup match area dmdtda 2000 to 2020:   {get_dmdtda(ds_dynamic_spinup_area):.2f}            kg m-2 yr-1')
print(f'Dynamical spinup match volume dmdtda 2000 to 2020: {get_dmdtda(ds_dynamic_spinup_volume):.2f}            kg m-2 yr-1')
print(f'Dynamical mu* calibration dmdtda 2000 to 2020:     {get_dmdtda(ds_dynamic_mu):.2f}            kg m-2 yr-1')

The blue curve shows the fixed geometry spinup (the old default spinup option). This option ignores area and length changes before the RGI date, which would be important, especially for hydrologic applications. Further, the short advance after the RGI date you see in the length is an “initial shock” after allowing the geometry to change. On the other side, the advantage of this method is that we match the area at the RGI date exact (red dot). And also the volume is matched perfectly because it is calculated backwards from the RGI date by using the constant area. The reference value of the length in these plots should not be considered as an observation as we are using an elevation band flowline, which loses some geometric information during the construction (see documentation).

Both dynamic spinup options (matching the area in green, matching the volume in red) address the issues of constant geometry and the “initial shock”. With this method, we get a changing surface geometry before the RGI date without an “initial Shock” afterwards. You also see clearly that depending on the variable to match, we are closer to the reference area or the reference volume. Which option to use depends on the use case and the available data. On a global scale, we selected to match to the area by default. A reason for this is that the area is a direct observation, whereas the volume is a consensus model derivation (Farinotti 2019).

Also shown in the plots above is the dynamic \(\mu^*\) calibration (including a dynamic spinup matching the area) in orange. However, you can not see an orange curve as it is exactly below the dynamic spinup matching the area (in green). The reason for this can be explained by comparing the reference geodetic mass balance (Hugonnet 2021) with the modelled ones, below the plots. You see that (by chance) the dynamic spinup already fits inside the provided observation uncertainty and therefore the dynamic \(\mu^*\) calibration does not need to look further. But to showcase the ability of the dynamic \(\mu^*\) calibration, we can artificially set the uncertainty lower and let it run again:

# define an artificial error for dmdtda
dmdtda_reference_error_artificial = 5 # error must be given as a positive number

tasks.run_dynamic_mu_star_calibration(gdir,
                                      ys=spinup_start_yr,  # When to start the spinup
                                      ye=2020,  # When the simulation should stop
                                      output_filesuffix='_dynamic_mu_artificial',  # Where to write the output
                                      ref_dmdtda=dmdtda_reference,  # user-provided geodetic mass balance observation
                                      err_ref_dmdtda=dmdtda_reference_error_artificial,  # uncertainty of user-provided geodetic mass balance observation 
                                     );

with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_dynamic_mu_artificial')) as ds:
    ds_dynamic_mu_artificial = ds.load()

# plot again
y0 = gdir.rgi_date + 1

f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 5))

ds_hist.volume_m3.plot(ax=ax1, label='Fixed geometry spinup');
ds_dynamic_mu_artificial.volume_m3.plot(ax=ax1, label='Dynamical mu* calibration');
ds_dynamic_spinup_area.volume_m3.plot(ax=ax1, label='Dynamical spinup match area');
ax1.set_title('Volume');
ax1.scatter(y0, volume_reference, c='C3', label='Reference values');
ax1.legend();

ds_hist.area_m2.plot(ax=ax2);
ds_dynamic_mu_artificial.area_m2.plot(ax=ax2);
ds_dynamic_spinup_area.area_m2.plot(ax=ax2);
ax2.set_title('Area');
ax2.scatter(y0, area_reference, c='C3')

ds_hist.length_m.plot(ax=ax3);
ds_dynamic_mu_artificial.length_m.plot(ax=ax3);
ds_dynamic_spinup_area.length_m.plot(ax=ax3);
ax3.set_title('Length');
ax3.scatter(y0, ds_hist.sel(time=y0).length_m, c='C3')

plt.tight_layout()
plt.show();

# and print out the modeled geodetic mass balances for comparision
def get_dmdtda(ds):
    yr0_ref_mb, yr1_ref_mb = ref_period.split('_')
    yr0_ref_mb = int(yr0_ref_mb.split('-')[0])
    yr1_ref_mb = int(yr1_ref_mb.split('-')[0])

    return ((ds.volume_m3.loc[yr1_ref_mb].values -
             ds.volume_m3.loc[yr0_ref_mb].values) /
            gdir.rgi_area_m2 /
            (yr1_ref_mb - yr0_ref_mb) *
            cfg.PARAMS['ice_density'])

print(f'Reference dmdtda 2000 to 2020 (artificial error):  {dmdtda_reference:.2f} +/- {dmdtda_reference_error_artificial:6.2f} kg m-2 yr-1')
print(f'Fixed geometry spinup dmdtda 2000 to 2020:         {get_dmdtda(ds_hist):.2f}            kg m-2 yr-1')
print(f'Dynamical spinup match area dmdtda 2000 to 2020:   {get_dmdtda(ds_dynamic_spinup_area):.2f}            kg m-2 yr-1')
print(f'Dynamical mu* calibration dmdtda 2000 to 2020:     {get_dmdtda(ds_dynamic_mu_artificial):.2f}            kg m-2 yr-1')

It is a made-up example, but you see here is a difference between the dynamic spinup and the dynamic \(\mu^*\) calibration. And indeed the dynamic mu star calibration manages to match our ‘very accurate’ geodetic mass balance.

For additional options of the two run tasks have a look at their doc-strings (run_dynamic_spinup and run_dynamic_mu_star_calibration). You can also change all the settings of the included dynamic spinup inside the run_dynamic_mu_star_calibration by using kwargs_run_function (e.g. kwargs_run_function = {'minimize_for':'volume'} to match for volume instead of area).

Example where it is not working so well#

The dynamic spinup however is not always working that well, as the example below shows:

# now we switch to the 'bad' glacier
gdir = gdirs[1]

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

# here we save the original mu_star for later
mu_star_original = gdir.read_json('local_mustar')['mu_star_glacierwide']

# and save some reference values we can compare to later
area_reference = gdir.rgi_area_m2  # RGI area
volume_reference = tasks.get_inversion_volume(gdir)  # volume after regional calibration to consensus
ref_period = cfg.PARAMS['geodetic_mb_period']  # get the reference period we want to use, 2000 to 2020
df_ref_dmdtda = utils.get_geodetic_mb_dataframe().loc[gdir.rgi_id]  # get the data from Hugonnet 2021
df_ref_dmdtda = df_ref_dmdtda.loc[df_ref_dmdtda['period'] == ref_period]  # only select the desired period
dmdtda_reference = df_ref_dmdtda['dmdtda'].values[0] * 1000  # get the reference dmdtda and convert into kg m-2 yr-1
dmdtda_reference_error = df_ref_dmdtda['err_dmdtda'].values[0] * 1000  # corresponding uncertainty

# ---- First do the fixed geometry spinup ----
tasks.run_from_climate_data(gdir,
                            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()

# ---- Second the dynamic spinup alone, matching area ----
tasks.run_dynamic_spinup(gdir,
                         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_area',  # Where to write the output
                         ye=2020,  # When the simulation should stop
                         );
# Read the output
with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_spinup_dynamic_area')) as ds:
    ds_dynamic_spinup_area = ds.load()
    
# ---- Third the dynamic spinup alone, matching volume ----
tasks.run_dynamic_spinup(gdir,
                         spinup_start_yr=spinup_start_yr,  # When to start the spinup
                         minimise_for='volume',  # what target to match at the RGI date
                         output_filesuffix='_spinup_dynamic_volume',  # Where to write the output
                         ye=2020,  # When the simulation should stop
                         );
# Read the output
with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_spinup_dynamic_volume')) as ds:
    ds_dynamic_spinup_volume = ds.load()

# Now make a plot for comparision
y0 = gdir.rgi_date + 1

f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 5))

ds_hist.volume_m3.plot(ax=ax1, label='Fixed geometry spinup');
ds_dynamic_spinup_area.volume_m3.plot(ax=ax1, label='Dynamical spinup match area');
ds_dynamic_spinup_volume.volume_m3.plot(ax=ax1, label='Dynamical spinup match volume');
ax1.set_title('Volume');
ax1.scatter(y0, volume_reference, c='C3', label='Reference values');
ax1.legend();

ds_hist.area_m2.plot(ax=ax2);
ds_dynamic_spinup_area.area_m2.plot(ax=ax2);
ds_dynamic_spinup_volume.area_m2.plot(ax=ax2);
ax2.set_title('Area');
ax2.scatter(y0, area_reference, c='C3')

ds_hist.length_m.plot(ax=ax3);
ds_dynamic_spinup_area.length_m.plot(ax=ax3);
ds_dynamic_spinup_volume.length_m.plot(ax=ax3);
ax3.set_title('Length');
ax3.scatter(y0, ds_hist.sel(time=y0).length_m, c='C3')

plt.tight_layout()
plt.show();

# and print out the modeled geodetic mass balances for comparision
def get_dmdtda(ds):
    yr0_ref_mb, yr1_ref_mb = ref_period.split('_')
    yr0_ref_mb = int(yr0_ref_mb.split('-')[0])
    yr1_ref_mb = int(yr1_ref_mb.split('-')[0])

    return ((ds.volume_m3.loc[yr1_ref_mb].values -
             ds.volume_m3.loc[yr0_ref_mb].values) /
            gdir.rgi_area_m2 /
            (yr1_ref_mb - yr0_ref_mb) *
            cfg.PARAMS['ice_density'])

print(f'Reference dmdtda 2000 to 2020 (Hugonnet 2021):     {dmdtda_reference:.2f} +/- {dmdtda_reference_error:6.2f} kg m-2 yr-1')
print(f'Fixed geometry spinup dmdtda 2000 to 2020:         {get_dmdtda(ds_hist):.2f}            kg m-2 yr-1')
print(f'Dynamical spinup match area dmdtda 2000 to 2020:   {get_dmdtda(ds_dynamic_spinup_area):.2f}            kg m-2 yr-1')
print(f'Dynamical spinup match volume dmdtda 2000 to 2020: {get_dmdtda(ds_dynamic_spinup_volume):.2f}            kg m-2 yr-1')

In this example, you see that the dynamic spinup runs do not start in 1980 but other years. The reason for this are the two main problems and the coping strategy of reducing the spinup time, described here in more detail. Be aware that even if you set a specific spinup year the resulting model run can start at a later year. If you want to get an error when this is happening you can use the min_spinup_period. For example you can set min_spinup_period = gdir.rgi_date - spinup_start_yr, with this an error is raised whenever it is not possible to run a dynamic spinup over the whole period.

Dynamical spinup and hydrological output#

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

# Here we use Khumbu again
gdir = gdirs[0]

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

# Dynamic mu star calibration including a dynamic spinup to match for area
tasks.run_with_hydro(gdir,
                     run_task=tasks.run_dynamic_mu_star_calibration,  # The task to run with hydro
                     ys=spinup_start_yr,  # When to start the spinup
                     ye=2020,  # When the simulation should stop
                     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
                     );

# Read the output
with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix='_spinup_dynamic_hydro')) as ds:
    ds_dynamic_spinup_hydro = ds.load()

ds_dynamic_spinup_hydro = ds_dynamic_spinup_hydro.isel(time=slice(0, -1))  # The last timestep is incomplete for hydro (not started)

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_dynamic_spinup_hydro.variables if 'month_2d' not in ds_dynamic_spinup_hydro[v].dims]
df_annual = ds_dynamic_spinup_hydro[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_ids[0]);

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'); plt.title(rgi_ids[0]);

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

What’s next?#

return to the OGGM documentation
back to the table of contents

OGGM tutorials

Dynamic spinup and dynamic \mu^* calibration for past simulations

Contents

Dynamic spinup and dynamic \(\mu^*\) calibration for past simulations#

High-level explanation of `run_dynamic_spinup`:#

Two main problems why the dynamic spinup could not work:#

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

High-level explanation of `run_dynamic_mu_star_calibration` :#

The minimization algorithm:#

Start with glacier example#

Define the glaciers we will play with#

Preparing the glacier data#

Start by comparing three different spinups#

Example where it is not working so well#

Dynamical spinup and hydrological output#

Annual runoff#

Monthly runoff#

What’s next?#

OGGM tutorials

Dynamic spinup and dynamic \mu^* calibration for past simulations

Contents

Dynamic spinup and dynamic \(\mu^*\) calibration for past simulations#

High-level explanation of run_dynamic_spinup:#

Two main problems why the dynamic spinup could not work:#

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

High-level explanation of run_dynamic_mu_star_calibration :#

The minimization algorithm:#

Start with glacier example#

Define the glaciers we will play with#

Preparing the glacier data#

Start by comparing three different spinups#

Example where it is not working so well#

Dynamical spinup and hydrological output#

Annual runoff#

Monthly runoff#

What’s next?#

High-level explanation of `run_dynamic_spinup`:#

High-level explanation of `run_dynamic_mu_star_calibration` :#