Ice thickness inversion

This example shows how to run the OGGM ice thickness inversion model with various ice parameters: the deformation parameter A and a sliding parameter (fs).

There is currently no “best” set of parameters for the ice thickness inversion model. As shown in Maussion et al. (2019), the default parameter set results in global volume estimates which are a bit larger than previous values. For the consensus estimate of Farinotti et al. (2019), OGGM participated with a deformation parameter A that is 1.5 times larger than the generally accepted default value.

There is no reason to think that the ice parameters are the same between neighboring glaciers. There is currently no “good” way to calibrate them, or at least no generaly accepted one. We won’t discuss the details here, but we provide a script to illustrate the sensitivity of the model to this choice.

New in version 1.4: we demonstrate how to apply a new global task in OGGM, workflow.calibrate_inversion_from_consensus to calibrate the A parameter to match the consensus estimate from Farinotti et al. (2019).

At the end of this tutorial, we show how to distribute the “flowline thicknesses” on a glacier map.

Run

# Libs
import geopandas as gpd

# Locals
import oggm.cfg as cfg
from oggm import utils, workflow, tasks, graphics

# Initialize OGGM and set up the default run parameters
cfg.initialize(logging_level='WARNING')
rgi_region = '11'  # Region Central Europe

# Local working directory (where OGGM will write its output)
WORKING_DIR = utils.gettempdir('OGGM_Inversion')
utils.mkdir(WORKING_DIR, reset=True)
cfg.PATHS['working_dir'] = WORKING_DIR

# This is useful here
cfg.PARAMS['use_multiprocessing'] = True

# RGI file
path = utils.get_rgi_region_file(rgi_region)
rgidf = gpd.read_file(path)

# Select the glaciers in the Pyrenees
rgidf = rgidf.loc[rgidf['O2Region'] == '2']

# Sort for more efficient parallel computing
rgidf = rgidf.sort_values('Area', ascending=False)

# Go - get the pre-processed glacier directories
# We start at level 3, because we need all data for the inversion
gdirs = workflow.init_glacier_directories(rgidf, from_prepro_level=3, prepro_border=10)

# Because of recent changes in OGGM not yet available in the preprocessed directories, 
# we re-run this task:
workflow.execute_entity_task(tasks.prepare_for_inversion, gdirs)

# Default parameters
# Deformation: from Cuffey and Patterson 2010
glen_a = 2.4e-24
# Sliding: from Oerlemans 1997
fs = 5.7e-20
2021-05-18 08:17:07: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2021-05-18 08:17:07: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2021-05-18 08:17:07: oggm.cfg: Multiprocessing: using all available processors (N=2)
2021-05-18 08:17:07: oggm.cfg: Multiprocessing switched ON after user settings.
2021-05-18 08:17:33: oggm.workflow: init_glacier_directories from prepro level 3 on 35 glaciers.
2021-05-18 08:17:33: oggm.workflow: Execute entity task gdir_from_prepro on 35 glaciers
2021-05-18 08:17:36: oggm.workflow: Execute entity task prepare_for_inversion on 35 glaciers
with utils.DisableLogger():  # this scraps some output - to use with caution!!!
    
    # Correction factors
    factors = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
    factors += [1.1, 1.2, 1.3, 1.5, 1.7, 2, 2.5, 3, 4, 5]
    factors += [6, 7, 8, 9, 10]

    # Run the inversions tasks with the given factors
    for f in factors:
        # Without sliding
        suf = '_{:03d}_without_fs'.format(int(f * 10))
        workflow.execute_entity_task(tasks.mass_conservation_inversion, gdirs,
                                     glen_a=glen_a*f, fs=0)
        # Store the results of the inversion only
        utils.compile_glacier_statistics(gdirs, filesuffix=suf,
                                         inversion_only=True)

        # With sliding
        suf = '_{:03d}_with_fs'.format(int(f * 10))
        workflow.execute_entity_task(tasks.mass_conservation_inversion, gdirs,
                                     glen_a=glen_a*f, fs=fs)
        # Store the results of the inversion only
        utils.compile_glacier_statistics(gdirs, filesuffix=suf,
                                         inversion_only=True)

Read the data

The data are stored as csv files in the working directory. The easiest way to read them is to use pandas!

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy import stats
import os

Let’s read the output of the inversion with the default OGGM parameters first:

df = pd.read_csv(os.path.join(WORKING_DIR, 'glacier_statistics_011_without_fs.csv'), index_col=0)
df
rgi_region rgi_subregion name cenlon cenlat rgi_area_km2 rgi_year glacier_type terminus_type is_tidewater ... inv_volume_km3 vas_volume_km3 inv_volume_bsl_km3 inv_volume_bwl_km3 dem_source n_orig_centerlines flowline_type perc_invalid_flowline inversion_glen_a inversion_fs
rgi_id
RGI60-11.03208 11 11-02 Aneto 0.646032 42.641357 0.622 2011 Glacier Land-terminating False ... 0.029854 0.017699 0.0 0.0 NASADEM 2 centerlines 0.205128 2.640000e-24 0
RGI60-11.03232 11 11-02 Ossoue -0.141222 42.771191 0.449 2011 Glacier Land-terminating False ... 0.019636 0.011306 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03222 11 11-02 Monte Perdido Lower 0.039655 42.679337 0.308 2011 Glacier Land-terminating False ... 0.009120 0.006733 0.0 0.0 NASADEM 2 centerlines 0.000000 2.640000e-24 0
RGI60-11.03209 11 11-02 Maladeta E 0.639343 42.648979 0.260 2011 Glacier Land-terminating False ... 0.007672 0.005334 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03231 11 11-02 Oulettes Du Gaube -0.144215 42.779716 0.130 2011 Glacier Land-terminating False ... 0.002747 0.002057 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03213 11 11-02 Seil De La Baque E 0.494304 42.694481 0.105 2011 Glacier Land-terminating False ... 0.002118 0.001533 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03239 11 11-02 Jezerces Gl No1 19.807627 42.465839 0.097 2007 Glacier Land-terminating False ... 0.002538 0.001375 0.0 0.0 NASADEM 1 centerlines 0.562500 2.640000e-24 0
RGI60-11.03229 11 11-02 Gabietous -0.057456 42.695129 0.087 2011 Glacier Land-terminating False ... 0.001720 0.001184 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03218 11 11-02 Llardana 0.428193 42.653946 0.084 2011 Glacier Land-terminating False ... 0.001472 0.001128 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03228 11 11-02 Taillon -0.039683 42.695419 0.083 2011 Glacier Land-terminating False ... 0.001733 0.001110 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03205 11 11-02 Tempestades 0.666092 42.626839 0.077 2011 Glacier Land-terminating False ... 0.001395 0.001001 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03206 11 11-02 Barrancs 0.658785 42.632011 0.067 2011 Glacier Land-terminating False ... 0.001078 0.000827 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03217 11 11-02 La Paul 0.437356 42.661072 0.057 2011 Glacier Land-terminating False ... 0.000993 0.000662 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03221 11 11-02 Monte Perdido Upper 0.033972 42.678589 0.054 2011 Glacier Land-terminating False ... 0.001139 0.000614 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03210 11 11-02 Maladeta W 0.632566 42.652107 0.054 2011 Glacier Land-terminating False ... 0.000884 0.000614 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03233 11 11-02 Infierno Central -0.258849 42.783802 0.053 2011 Glacier Land-terminating False ... 0.000868 0.000599 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03225 11 11-02 Astazou 0.034118 42.703136 0.052 2011 Glacier Land-terminating False ... 0.001376 0.000583 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03211 11 11-02 Boum 0.555278 42.700157 0.050 2011 Glacier Land-terminating False ... 0.001101 0.000553 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03240 11 11-02 Jezerces Gl No2 19.815450 42.447353 0.049 2007 Glacier Land-terminating False ... 0.000666 0.000538 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03227 11 11-02 Pailla W 0.015287 42.705482 0.046 2011 Glacier Land-terminating False ... 0.000711 0.000493 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03219 11 11-02 Barroude 0.142224 42.726147 0.042 2011 Glacier Land-terminating False ... 0.000828 0.000435 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03220 11 11-02 Munia 0.129816 42.717735 0.041 2011 Glacier Land-terminating False ... 0.000583 0.000421 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03212 11 11-02 Portillon D Oo 0.511460 42.693188 0.040 2011 Glacier Land-terminating False ... 0.000549 0.000407 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03226 11 11-02 Pailla E 0.021260 42.705322 0.035 2011 Glacier Land-terminating False ... 0.000639 0.000339 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03241 11 11-02 It4L30000001 Gh Del Calderone 13.566710 42.470540 0.033 1990 Glacier Land-terminating False ... 0.000619 0.000312 0.0 0.0 NASADEM 1 centerlines 0.250000 2.640000e-24 0
RGI60-11.03230 11 11-02 Petit Vignemale -0.137718 42.775856 0.027 2011 Glacier Land-terminating False ... 0.000370 0.000237 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03215 11 11-02 Gourgs Blancs 0.480293 42.700859 0.026 2011 Glacier Land-terminating False ... 0.000468 0.000225 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03234 11 11-02 Las Neous -0.286945 42.838642 0.025 2011 Glacier Land-terminating False ... 0.000465 0.000213 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03223 11 11-02 Marbore 0.018435 42.693417 0.023 2011 Glacier Land-terminating False ... 0.000325 0.000190 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03224 11 11-02 Cilindro 0.026603 42.689346 0.022 2011 Glacier Land-terminating False ... 0.000319 0.000179 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03204 11 11-02 Valier 1.088501 42.798626 0.021 2011 Glacier Land-terminating False ... 0.000264 0.000168 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03214 11 11-02 Seil De La Baque W 0.490498 42.696220 0.020 2011 Glacier Land-terminating False ... 0.000303 0.000157 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03238 11 11-02 Debeli Namet Gl 19.067500 43.115000 0.019 2003 Glacier Land-terminating False ... 0.000180 0.000146 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03207 11 11-02 Coronas 0.654112 42.632977 0.015 2011 Glacier Land-terminating False ... 0.000259 0.000106 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0
RGI60-11.03216 11 11-02 Posets 0.438703 42.656570 0.010 2011 Glacier Land-terminating False ... 0.000130 0.000060 0.0 0.0 NASADEM 1 centerlines 0.000000 2.640000e-24 0

35 rows × 21 columns

There are only 35 glaciers in the Pyrenees! That’s why the run was relatively fast.

Results

One way to visualize the output is to plot the volume as a function of area in a log-log plot, illustrating the well known volume-area relationship of mountain glaciers:

ax = df.plot(kind='scatter', x='rgi_area_km2', y='inv_volume_km3')
ax.semilogx(); ax.semilogy()
xlim, ylim = [1e-2, 0.7], [1e-5, 0.05]
ax.set_xlim(xlim); ax.set_ylim(ylim);
../_images/inversion_13_0.png

As we can see, there is a clear relationship, but it is not perfect. Let’s fit a line to these data (the “volume-area scaling law”):

# Fit in log space 
dfl = np.log(df[['inv_volume_km3', 'rgi_area_km2']])
slope, intercept, r_value, p_value, std_err = stats.linregress(dfl.rgi_area_km2.values, dfl.inv_volume_km3.values)

In their seminal paper, Bahr et al. (1997) describe this relationship as:

\[V = \alpha S^{\gamma}\]

With V the volume in km\(^3\), S the area in km\(^2\) and \(\alpha\) and \(\gamma\) the scaling parameters (0.034 and 1.375, respectively). How does OGGM compare to these in the Pyrenees?

print('power: {:.3f}'.format(slope))
print('slope: {:.3f}'.format(np.exp(intercept)))
power: 1.301
slope: 0.045
ax = df.plot(kind='scatter', x='rgi_area_km2', y='inv_volume_km3', label='OGGM glaciers')
ax.plot(xlim, np.exp(intercept) * (xlim ** slope), color='C3', label='Fitted line')
ax.semilogx(); ax.semilogy()
ax.set_xlim(xlim); ax.set_ylim(ylim);
ax.legend();
../_images/inversion_18_0.png

Sensitivity analysis

Now, let’s read the output files of each run separately, and compute the regional volume out of them:

dftot = pd.DataFrame(index=factors)
for f in factors:
    # Without sliding
    suf = '_{:03d}_without_fs'.format(int(f * 10))
    fpath = os.path.join(WORKING_DIR, 'glacier_statistics{}.csv'.format(suf))
    _df = pd.read_csv(fpath, index_col=0, low_memory=False)
    dftot.loc[f, 'without_sliding'] = _df.inv_volume_km3.sum()
    
    # With sliding
    suf = '_{:03d}_with_fs'.format(int(f * 10))
    fpath = os.path.join(WORKING_DIR, 'glacier_statistics{}.csv'.format(suf))
    _df = pd.read_csv(fpath, index_col=0, low_memory=False)
    dftot.loc[f, 'with_sliding'] = _df.inv_volume_km3.sum()

And plot them:

dftot.plot();
plt.xlabel('Factor of Glen A (default 1)'); plt.ylabel('Regional volume (km$^3$)');
../_images/inversion_23_0.png

As you can see, there is quite a difference between the solutions. In particular, close to the default value for Glen A, the regional estimates are very sensitive to small changes in A. The calibration of A is a problem that has yet to be resolved by global glacier models…

New in version 1.4: calibrate to match the consensus estimate

Here, a “best Glen A” is found in order that the total inverted volume of the glaciers of gdirs fits to the 2019 consensus estimate.

cdf = workflow.calibrate_inversion_from_consensus(gdirs, filter_inversion_output=False)
2021-05-18 08:20:31: oggm.workflow: Applying global task calibrate_inversion_from_consensus on 35 glaciers
2021-05-18 08:20:31: oggm.workflow: Consensus estimate optimisation with A factor: 0.1 and fs: 0
2021-05-18 08:20:31: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2021-05-18 08:20:31: oggm.workflow: Execute entity task prepare_for_inversion on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Execute entity task get_inversion_volume on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Consensus estimate optimisation with A factor: 10.0 and fs: 0
2021-05-18 08:20:32: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Execute entity task prepare_for_inversion on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Execute entity task get_inversion_volume on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Consensus estimate optimisation with A factor: 9.927344372640782 and fs: 0
2021-05-18 08:20:32: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2021-05-18 08:20:32: oggm.workflow: Execute entity task prepare_for_inversion on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Execute entity task get_inversion_volume on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Consensus estimate optimisation with A factor: 9.503932371199765 and fs: 0
2021-05-18 08:20:33: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Execute entity task prepare_for_inversion on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Execute entity task get_inversion_volume on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Consensus estimate optimisation with A factor: 9.551452033056764 and fs: 0
2021-05-18 08:20:33: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2021-05-18 08:20:33: oggm.workflow: Execute entity task prepare_for_inversion on 35 glaciers
2021-05-18 08:20:34: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers
2021-05-18 08:20:34: oggm.workflow: Execute entity task get_inversion_volume on 35 glaciers
2021-05-18 08:20:34: oggm.workflow: calibrate_inversion_from_consensus converged after 4 iterations and fs=0. The resulting Glen A factor is 9.503932371199765.
2021-05-18 08:20:34: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2021-05-18 08:20:34: oggm.workflow: Execute entity task prepare_for_inversion on 35 glaciers
2021-05-18 08:20:34: oggm.workflow: Execute entity task mass_conservation_inversion on 35 glaciers
2021-05-18 08:20:34: oggm.workflow: Execute entity task get_inversion_volume on 35 glaciers
cdf.sum()
vol_itmix_m3        6.423913e+07
vol_bsl_itmix_m3    0.000000e+00
vol_oggm_m3         6.425211e+07
dtype: float64

Note that here we calibrate the Glen A parameter to a value that is equal for all glaciers of gdirs (here \(A \sim 9.504\cdot A_0\)), i.e. we calibrate to match the total volume of all glaciers and not to match them individually.

cdf.iloc[:3]
vol_itmix_m3 vol_bsl_itmix_m3 vol_oggm_m3
RGIId
RGI60-11.03208 1.564802e+07 0.0 1.972692e+07
RGI60-11.03232 1.422478e+07 0.0 1.340392e+07
RGI60-11.03222 5.158277e+06 0.0 6.604405e+06

just as a side note, “vol_bsl_itmix_m3” means volume below sea level and is therefore zero for these Alpine glaciers!

Distributed ice thickness

The OGGM inversion and dynamical models use the “1D” flowline assumption: for some applications, you might want to use OGGM to create distributed ice thickness maps. Currently, OGGM implements two ways to “distribute” the flowline thicknesses, but only the simplest one works robustly:

# Distribute
workflow.execute_entity_task(tasks.distribute_thickness_per_altitude, gdirs);
2021-05-18 08:20:34: oggm.workflow: Execute entity task distribute_thickness_per_altitude on 35 glaciers

We just created a new output of the model, which we can access in the gridded_data file:

# xarray is an awesome library! Did you know about it?
import xarray as xr
ds = xr.open_dataset(gdirs[0].get_filepath('gridded_data'))
ds.distributed_thickness.plot();
../_images/inversion_36_0.png

Plot many glaciers on a map

Let’s select a group of glaciers close to each other:

rgi_ids = ['RGI60-11.0{}'.format(i) for i in range(3205, 3211)]
sel_gdirs = [gdir for gdir in gdirs if gdir.rgi_id in rgi_ids]
graphics.plot_googlemap(sel_gdirs)
# you might need to install motionless if it is not yet in your environment
../_images/inversion_39_0.png

Using OGGM

Since a recent PR (21.05.2020), OGGM can plot the thickness of a group of glaciers on a map:

graphics.plot_distributed_thickness(sel_gdirs)
../_images/inversion_42_0.png

This is however not very useful because OGGM can only plot on a map as large as the local glacier map of the first glacier in the list. See this issue for a discussion about why.

Using salem

Under the hood, OGGM uses salem to make the plots. Let’s do that for our case: it requires some manual tweaking, but it should be possible to automatize this better in the future.

Note: this also requires a very recent version of salem to work (21.05.2020)

import salem
# Make a grid covering the desired map extent
g = salem.mercator_grid(center_ll=(0.65, 42.64), extent=(4000, 4000))
# Create a map out of it
smap = salem.Map(g, countries=False)
# Add the glaciers outlines
for gdir in sel_gdirs:
    crs = gdir.grid.center_grid
    geom = gdir.read_pickle('geometries')
    poly_pix = geom['polygon_pix']
    smap.set_geometry(poly_pix, crs=crs, fc='none', zorder=2, linewidth=.2)
    for l in poly_pix.interiors:
        smap.set_geometry(l, crs=crs, color='black', linewidth=0.5)
f, ax = plt.subplots(figsize=(6, 6))
smap.visualize();
../_images/inversion_46_0.png
# Now add the thickness data
for gdir in sel_gdirs:
    grids_file = gdir.get_filepath('gridded_data')
    with utils.ncDataset(grids_file) as nc:
        vn = 'distributed_thickness'
        thick = nc.variables[vn][:]
        mask = nc.variables['glacier_mask'][:]
    thick = np.where(mask, thick, np.NaN)
    # The "overplot=True" is key here
    # this needs a recent version of salem to run properly
    smap.set_data(thick, crs=gdir.grid.center_grid, overplot=True)
# Set colorscale and other things
smap.set_cmap(graphics.OGGM_CMAPS['glacier_thickness'])
smap.set_plot_params(nlevels=256)
# Plot
f, ax = plt.subplots(figsize=(6, 6))
smap.visualize(ax=ax, cbar_title='Glacier thickness (m)');
../_images/inversion_48_0.png

What’s next?