Understand the difference between the ice dynamic solvers in OGGM

In version 1.6, OGGM changed the default numeric solver to the Semi-Implicit model. In this notebook, we explore the main differences compared to the old default, the Flux-Based model.

import time
import xarray as xr
import matplotlib.pyplot as plt
from oggm import cfg, utils, workflow, graphics, tasks
from oggm.core.flowline import FluxBasedModel, SemiImplicitModel

# Initialize OGGM and set up the default run parameters
cfg.initialize(logging_level='WARNING')

# Define our test glacier (Baltoro)
rgi_ids = ['RGI60-14.06794']

# load elevation band representation
cfg.PATHS['working_dir'] = utils.gettempdir('OGGM_dynamic_solvers_elevation_bands', reset=True)
base_url_eb = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/elev_bands/W5E5/'
gdir_eb = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_base_url=base_url_eb)[0]

# load centerline representation
cfg.PATHS['working_dir'] = utils.gettempdir('OGGM_dynamic_solvers_centerliens', reset=True)
base_url_cl = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/centerlines/W5E5/'
gdir_cl = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_base_url=base_url_cl)[0]

2024-08-25 21:28:16: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-08-25 21:28:16: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-08-25 21:28:16: oggm.cfg: Multiprocessing: using all available processors (N=4)
2024-08-25 21:28:17: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2024-08-25 21:28:17: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
2024-08-25 21:28:18: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2024-08-25 21:28:18: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

Flux-Based model is more flexible, but unstable

The big advantage of the Flux-Based model is that it works for all flowline representations (multiple flowlines and different bed shapes). See the tutorial “elevation band” and “centerline” flowlines for a hands-on introduction to the different flowline types.

# run Flux-Based with centerlines
tasks.run_random_climate(gdir_cl,
                         evolution_model=FluxBasedModel,
                         nyears=300,
                         y0=2000,
                         seed=0,
                         store_fl_diagnostics=True,
                         output_filesuffix='_flux_based')

# plot result
with xr.open_dataset(gdir_cl.get_filepath('model_diagnostics', filesuffix='_flux_based')) as ds:
    ds_trap = ds.load()
ds_trap.volume_m3.plot();

# run Flux-Based with elevation bands
start_time = time.time()  # time it for later comparision
tasks.run_random_climate(gdir_eb,
                         evolution_model=FluxBasedModel,
                         nyears=300,
                         y0=2000,
                         seed=0,
                         store_fl_diagnostics=True,
                         output_filesuffix='_flux_based')
flux_based_time = time.time() - start_time

# plot result
with xr.open_dataset(gdir_eb.get_filepath('model_diagnostics', filesuffix='_flux_based')) as ds:
    ds_flux_eb = ds.load()
ds_flux_eb.volume_m3.plot();

Whereas the Semi-Impicit model only works for single trapezoidal flowlines (elevation bands).

# run Semi-Implicit with centerlines raises an error 
# The code below would fail (expected)
import pytest
with pytest.raises(ValueError):
    tasks.run_random_climate(gdir_cl,
                             evolution_model=SemiImplicitModel,
                             y0=2000,
                             seed=0,
                             store_fl_diagnostics=True,
                             output_filesuffix='_semi_implicit')

2024-08-25 21:31:16: oggm.core.flowline: ValueError occurred during task flowline_model_run_semi_implicit on RGI60-14.06794: Implicit model does not work with tributaries.
2024-08-25 21:31:16: oggm.core.flowline: ValueError occurred during task run_random_climate_semi_implicit on RGI60-14.06794: Implicit model does not work with tributaries.

# run Semi-Implicit with elevation bands
start_time = time.time()  # time it for later comparision
tasks.run_random_climate(gdir_eb,
                         evolution_model=SemiImplicitModel,
                         nyears=300,
                         y0=2000,
                         seed=0,
                         store_fl_diagnostics=True,
                         output_filesuffix='_semi_implicit')
semi_implicit_time = time.time() - start_time

# plot result
with xr.open_dataset(gdir_eb.get_filepath('model_diagnostics', filesuffix='_semi_implicit')) as ds:
    ds_impl_eb = ds.load()

ds_impl_eb.volume_m3.plot(label='SemiImplicitModel', lw=4)
ds_flux_eb.volume_m3.plot(label='FluxBasedModel')
plt.legend();

You see that for the elevation band flowlines, both produce similar results. The differences arise from numeric instabilities in the Flux-Based model (see next paragraph). You can redo the experiment with a glacier where these instabilities are not that severe (e.g. RGI60-11.00897 Hintereisferner) and you will see both models produce the same result.

Semi-Implicit model is faster and more stable, but less flexible

Even the Semi-Implicit model is not as flexible as the Flux-Based one, we see it is faster when comparing the computing time:

print(f'Semi-Implicit time needed: {semi_implicit_time:.1f} s')
print(f'Flux-Based time needed: {flux_based_time:.1f} s')

Semi-Implicit time needed: 1.8 s
Flux-Based time needed: 3.0 s

For a single glacier, this speed-up is probably not that important, but when thinking about regional to global simulations it can save you a lot of time.

One reason for the speed-up is that the Semi-Implicit model is numerically more stable and can take larger time steps without producing instabilities:

# open flowline diagnostics
f_impl = gdir_eb.get_filepath('fl_diagnostics', filesuffix='_semi_implicit')
f_flux = gdir_eb.get_filepath('fl_diagnostics', filesuffix='_flux_based')
with xr.open_dataset(f_impl, group=f'fl_0') as ds:
    ds_fl_impl = ds.load()
with xr.open_dataset(f_flux, group=f'fl_0') as ds:
    ds_fl_flux = ds.load()
    
# compare velocities along flowline
year = 100
ds_fl_impl.sel(time=year).ice_velocity_myr.plot(label='SemiImplicitModel')
ds_fl_flux.sel(time=year).ice_velocity_myr.plot(label='FluxBasedModel')
plt.legend();

In this case instabilities are visible for the FluxBasedModel at around 30 km distance along the flowline. They can lead to very large velocities which reduce the maximum possible step size due to the cfl-criterion (see also in the documentation).

The increased computational speed and, even more importantly, the increased stability are the reasons why we switched to the SemiImplicitModel in OGGM v1.6.

However, if you want to set the FluxBasedModel as your default, you can do so with:

cfg.PARAMS['evolution_model'] = 'FluxBased'  # default is 'SemiImplicit'

2024-08-25 21:31:18: oggm.cfg: PARAMS['evolution_model'] changed from `SemiImplicit` to `FluxBased`.

Have 5 minutes more? The bed shape of the downstream line

This paragraph deals with the downstream line, the initially ice-free part in front of the glacier. You can see it below as the red line connecting the end of the outline with the left border of the figure:

graphics.plot_centerlines(gdir_cl,
                          use_flowlines=True,
                          add_downstream=True)

In OGGM before v1.6, with the FluxBasedModel, the shape of this downstream line was defined by fitting a parabola to the valley walls. However, for the SemiImplicitModel we had to change the shape to a trapezoidal, eventhough a parabola approximates a mountain valley arguably better. We checked the influence of this change on advancing glaciers and found negligibly small differences in the volume on a regional scale. There might be some differences in the area.

By default, we use a trapezoidal bed shape for the downstream line:

fl_trap = gdir_eb.read_pickle('model_flowlines')
fl_trap[-1].is_trapezoid[fl_trap[-1].thick == 0]

array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True])

But if for any reason you decided to use the FluxBasedModel you also can switch back to a parabolic downstream line using cfg.PARAMS['downstream_line_shape'] = 'parabola'.

# change the downstream line shape
cfg.PARAMS['downstream_line_shape'] = 'parabola'  # default is 'trapezoidal'

# IMPORTANT: need to call init_present_time_glacier to take effect
tasks.init_present_time_glacier(gdir_eb)

fl_trap = gdir_eb.read_pickle('model_flowlines')
fl_trap[-1].is_trapezoid[fl_trap[-1].thick == 0]

2024-08-25 21:31:23: oggm.cfg: PARAMS['downstream_line_shape'] changed from `trapezoidal` to `parabola`.

array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False])

What’s next?

• return to the OGGM documentation

• back to the table of contents