v0.1 (29 March 2016)¶

Initial release, used to prepare the data submitted to ITMIX (see ITMIX Experiment 2016).

This release is the result of several months of development (outside of GitHub for a large part). Several people have contributed to this release:

• Michael Adamer (intern, UIBK), participated to the development of the centerline determination algorithm (2014)
• Kévin Fourteau (intern, UIBK, ENS Cachan), participated to the development of the inversion and the flowline modelling algorithms (2014-2015)
• Alexander H. Jarosch (Associate Professor, University of Iceland), developed the MUSCL-SuperBee model (PR23)
• Johannes Landmann (intern, UIBK), participated to the links between databases project (2015)
• Ben Marzeion (project leader, University of Bremen), scientific advisor
• Fabien Maussion (project leader, UIBK), core developer
• Felix Oesterle (Post-Doc, UIBK), develops OGGR and provided the AWS deployment script (PR25)
• Timo Rothenpieler (programmer, University of Bremen), participated to the OGGM deployment script (e.g. PR34, PR48), and developed OGGM installation tools
• Christian Wild (master student, UIBK), participated to the development of the centerline determination algorithm (2014)

Dependencies¶

Standard SciPy track:

• numpy
• scipy
• scikit-image
• pillow
• matplotlib
• pandas
• xarray
• joblib

Python 2 support:

• six

Configuration file parsing tool:

• configobj

I/O:

• netcdf4

GIS tools:

• gdal
• shapely
• pyproj
• rasterio
• geopandas

Testing:

• nose

Other libraries:

• salem
• cleo
• motionless (py3)

With conda (all platforms)¶

$ sudo apt-get install build-essential liblapack-dev gfortran libproj-dev gdal-bin libgdal-dev netcdf-bin ncview python-netcdf ttf-bitstream-vera

conda create --name oggm python=3.4

source activate oggm

conda install -c ioos geopandas matplotlib Pillow joblib netCDF4 scikit-image configobj nose pyproj numpy krb5 rasterio xarray

pip install git+https://github.com/fmaussion/motionless.git
pip install git+https://github.com/fmaussion/salem.git
pip install git+https://github.com/fmaussion/cleo.git

git clone git@github.com:OGGM/oggm.git

cd oggm

pip install -e .

nosetests .

With virtualenv (linux/debian)¶

For Debian / Ubuntu / Mint users only!

$ sudo apt-get install build-essential python-pip liblapack-dev gfortran libproj-dev

$ sudo apt-get install python-gtk2-dev

$ sudo apt-get install tk-dev python3-tk python3-dev

$ sudo apt-get install gdal-bin libgdal-dev python-gdal

$ sudo apt-get install netcdf-bin ncview python-netcdf

$ sudo pip install virtualenvwrapper

$ mkdir ~/.pyvirtualenvs

# Virtual environment options
export WORKON_HOME=$HOME/.pyvirtualenvs
source /usr/local/bin/virtualenvwrapper_lazy.sh

$ . ~/.profile

$ mkvirtualenv oggm_env -p /usr/bin/python

$ mkvirtualenv oggm_env -p /usr/bin/python3

$ workon oggm_env

$ pip install numpy scipy pandas shapely

$ ln -sf /usr/lib/python2.7/dist-packages/{glib,gobject,cairo,gtk-2.0,pygtk.py,pygtk.pth} $VIRTUAL_ENV/lib/python2.7/site-packages
$ pip install matplotlib

$ pip install matplotlib

>>> import matplotlib.pyplot as plt
>>> plt.plot([1,2,3])
>>> plt.show()

$ dpkg -s libgdal-dev

$ pip install gdal==1.10.0 --install-option="build_ext" --install-option="--include-dirs=/usr/include/gdal"
$ pip install fiona --install-option="build_ext" --install-option="--include-dirs=/usr/include/gdal"

$ pip install pyproj rasterio Pillow geopandas netcdf4 scikit-image configobj joblib xarray

$ pip install git+https://github.com/fmaussion/motionless.git
$ pip install git+https://github.com/fmaussion/salem.git
$ pip install git+https://github.com/fmaussion/cleo.git

Basics¶

Let \(S\) be the area of a cross-section perpendicular to the flowline. It has a width \(w\) and a thickness \(h\) and, in this example, a parabolic bed shape.

Example of a cross-section along the glacier flowline. Background image from http://www.swisseduc.ch/glaciers/alps/hintereisferner/index-de.html

Volume conservation for this discrete element implies:

where \(\dot{m}\) is the mass-balance, \(q = u S\) the flux of ice, and \(u\) the depth-integrated ice velocity ([Cuffey_Paterson_2010], p 310). This velocity can be computed from Glen’s flow law as a function of the basal shear stress \(\tau\):

The second term is to account for basal sliding, see e.g. [Oerlemans_1997] or [Golledge_Levy_2011]. It introduces an additional free parameter \(f_s\) and will therefore be ignored in a first approach. The deformation parameter \(f_d\) is better constrained and relates to Glen’s temperature‐dependent creep parameter \(A\):

The basal shear stress \(\tau\) depends e.g. on the geometry of the bed [Cuffey_Paterson_2010]. Currently it is assumed to be equal to the driving stress \(\tau_d\):

where \(\alpha\) is the slope of the flowline and \(\rho\) the density of ice. Both the FluxBasedModel and the MUSCLSuperBeeModel solve for these equations, but with different numerical schemes.

Bed thickness inversion¶

To compute the initial ice thikness \(h_0\), OGGM follows a methodology largely inspired from [Farinotti_etal_2009] but using a different apparent mass-balance (see also: Mass-balance) and another calibration algorithm.

The principle is simple. Let’s assume for now that we know the ice velocity \(u\) along the flowline of our present-time glacier. Then the above equations can be used to compute the section area \(S\) out of \(u\) and the other ice-flow parameters. Since we know the present-time width \(w\) with accuracy, \(h_0\) can be obtained by assuming a certain geometrical shape for the bed.

In OGGM, a number of climate and glacier related parameters are fixed prior to the inversion, leaving only one free parameter for the calibration of the bed inversion procedure: the inversion factor \(f_{inv}\). It is defined such as:

With \(A_0\) the standard creep parameter (2.4e-24). Currently, \(f_{inv}\) is calibrated to minimize the volume RMSD of all glaciers with a volume estimation in the GlaThiDa database. It is therefore neither glacier nor temperature dependent and does not account for uncertainties in GlaThiDa’s glacier-wide thickness estimations, two approximations which should be better handled in the future.

Flux based model¶

Most flowline models treat the volume conservation equation as a diffusion problem, taking advantage of the robust numerical solutions available for this type of equations. The problem with this approach is that it develops the \(\partial S / \partial t\) term to solve for ice thikness \(h\) directly, thus implying different diffusion equations for different bed geometries (e.g. [Oerlemans_1997] with a trapezoidal bed).

The OGGM flux based model solves for the \(\nabla \cdot q\) term on a staggered grid (hence the name). It has the advantage that the model numerics are the same for any bed shape, but it makes one important simplification: the stress \(\tau = \alpha \rho g h\) is always the same, regardless of the bed shape. This doesn’t mean that the shape has no influence on the glacier evolution, as explained below.

Glacier bed shapes¶

OGGM implements a number of possible bed-shapes. Currently the shape has no direct influence on ice dynamics, but it does influence how the width of the glacier changes with ice thickness and thus will influence the mass-balance \(w \, \dot{m}\). It appears that the flowline model is quite sensitive to the bed shape.

MUSCLSuperBeeModel¶

A shallow ice model with improved numerics ensuring mass-conservation in very steep walls [Jarosch_etal_2013]. The model is currently in development to account for various bed shapes and tributaries and will likely become the default in OGGM.

Glacier tributaries¶

Glaciers in OGGM have a main centerline and, sometimes, one or more tributaries (which can themsleves also have tributaries, see Centerlines). The number of these tributaries depends on many factors, but most of the time the algorithm works well.

The main flowline and its tributaries are all handled the same way and are modelled individually. The difference is that tributaries can transport mass to the branch they are flowing to. Numerically, this mass transport is handled by adding an element at the end of the flowline with the same properties (with, thickness...) as the last grid point, with the difference that the slope \(\alpha\) is computed with respect to the altitude of the point they are flowing to. The ice flux is then computed normally and transferred to the downlying branch.

The computation of the ice flux is always done first from the lowest order branches (without tributaries) to the highest ones, ensuring a correct mass-redistribution. The angle between tributary and main branch ensures that the former is not decoupled from the latter. If the angle is positive or if no ice is present at the end of the tributary, no mass exchange occurs.

References¶

Glacier divides¶

The glacier outlines provided by the RGI are the result of an automated (sometimes manually corrected) delineation procedure. In some cases, the RGI glacier is not the outline of a single glacier (in the dynamical sense) but is the outline of a “glacier complex”. OGGM does not (yet) implement a method to divide these glacier complexes into individual glaciers [1], but it does have a framework to deal with user provided glacier divides.

An RGI entity can have one or more “divides”, stored as polygons in a separate shapefile (oggm.cfg.set_divides_db()). Several OGGM tasks (such as compute_centerlines()) can work properly on single glaciers only, as shown on the example below (left: with divides, right: without).

The divides are handled as sub-directories in the glacier directory: each sub-directory stores the files specific to that divide. All the glaciers have at least one divide (“divide_01”), the divide ID 0 being reserved for the root directory which contains files concerning the entire RGI glacier. The tasks that work on divides only can be recognized by their div_id=None keyword argument (they also implement the divide_task decorator).

cfg.BASENAMES¶

This is a list of the files that can be found in the glacier directory or its divides:

apparent_mb.nc

The apparent mass-balance data needed for the inversion.

catchment_indices.pkl

A list of len n_centerlines, each element conaining a numpy array of the indices in the glacier grid which represent the centerline’s catchment area.

centerlines.pkl

A list of :py:class:oggm.Centerline instances, sorted by flow order.

climate_monthly.nc

The monthly climate timeseries for this glacier, stored in a netCDF file.

dem.tif

A geotiff file containing the DEM (reprojected into the local grid).

dem_source.pkl

A string with the source of the topo file (ASTER, SRTM, ...).

downstream_line.pkl

A shapely.LineString of the coordinates of the downstream line (flowing out of the glacier until the border of the domain) for each divide.

find_initial_glacier_params.pkl

geometries.pkl

A dict containing the shapely.Polygons of a divide. The “polygon_hr” entry contains the geometry transformed to the local grid in (i, j) coordinates, while the “polygon_pix” entry contains the geometries transformed into the coarse grid (the i, j elements are integers). The “polygon_area” entry contains the area of the polygon as computed by Shapely (it is needed because the divides will have their own area which is not obtained from the RGI file).

glacier_grid.pkl

A salem.Grid handling the georeferencing of the local grid.

gridded_data.nc

A netcdf file containing several gridded data variables such as topography, the glacier masks and more (see the netCDF file metadata).

inversion_flowlines.pkl

A “better” version of the Centerlines, now on a regular spacing i.e., not on the gridded (i, j) indices. The tails of the tributaries are cut out to make more realistic junctions. They are now “1.5D” i.e., with a width.

inversion_input.pkl

List of dicts containing the data needed for the inversion.

inversion_output.pkl

List of dicts containing the output data from the inversion.

inversion_params.pkl

Dict of fs and fd as computed by the inversion optimisation.

local_mustar.csv

A csv with three values: the local scalars mu*, t*, bias

major_divide.pkl

A simple integer in the glacier root directory (divide 00) containing the ID of the “major divide”, i.e. the one really flowing out of the glacier (the other downstream lines flowing into the main branch).

model_flowlines.pkl

List of flowlines ready to be run by the model.

mu_candidates.pkl

A pandas.Series with the (year, mu) data.

outlines.shp

The glacier outlines in the local projection.

past_model.pkl

past_models.pkl

ITMIX preprocessing¶

The glaciers are heterogeneous: valley glaciers, ice-caps, marine-terminating...

Blue: ITMIX outlines, Black: RGI outlines. White means that ITMIX didn’t provide the topography.

The first step that we needed to do is to formalize all this data so that OGGM can deal with it. This turned out to be a bit complicated since the ITMIX data was not standardized:

• we updated the RGI outlines with the ITMIX ones where possible
• for ice-caps, we kept the RGI outlines because OGGM currently needs the “pieces of cake” to compute the centerlines - see the important implications below.
• we decided to to compute the inversion on the OGGM local grids, not on the IMIX maps. This is important because OGGM makes decisions about the grid spacing it uses. Furthermore, the entire workflow is depending on these standardized local maps. We introduced a new routine in the pipeline in order to update the SRTM/ASTER topography with the ITMIX one. This also was not trivial because some ITMIX topographies are stopping directly at the glacier boundary.

OGGM preprocessing¶

Topography¶

OGGM uses following DEM data:

• SRTM V4.1 for [60S, 60N] (http://srtm.csi.cgiar.org/)
• GIMP DEM for Greenland (https://bpcrc.osu.edu/gdg/data/gimpdem)
• Corrected DEMs for Svalbard, Iceland, and North Canada (http://viewfinderpanoramas.org/dem3.html)
• ASTER V2 elsewhere

The corrected DEMs where necessary because ASTER data has many issues over glaciers. Take for example the DEM for two glaciers in Iceland:

Note that the hypsometry provided in RGI V5 also contains these errors. While the problems with the right plot are obvious, the glacier on the left (Dyngjujoekull) is practically impossible to filter automatically. On the plot below, I show the hypsometry that OGGM computed and the one by Mathias Huss:

Up to a few discrepancies due to projection issues, we both have the problem of non-zero bins below 750 m a.s.l. Fortunately, thanks to the work by Jonathan de Ferranti, these problems are now resolved in OGGM:

There is potential for even better coverage of corrected DEM, but this would require a bit more work (J. de Ferranti’s data is not always logically structured).

Calibration data¶

Another obstacle to global coverage was that the databases required for calibration (WGMS FoG and GlaThiDa) are not “linked” to RGI, i.e. there is no way to know which RGI entity corresponds to each database entry. Thanks to the work of Johannes Landmann at UIBK we now have comprehensive links with global coverage.

Some of the ITMIX glaciers were in the database, I removed them manually for the sake of the experiment (“blind run”):

• WGMS: Kesselwandferner, Brewster, Devon, Elbrus, Freya, Hellstugubreen, Urumqi
• GlaThiDa: Kesselwandferner, Unteraar

These leaves us with 201861 WGMS glaciers with at least 5 years of mass-balance data available for calibration of the mass-balance:

And 133 GlaThiDa glaciers with glacier-wide average thickness estimates. The coverage of GlaThiDa is not very good, which is probably a problem:

Inversion procedure¶

Refer to the general documentation for details about the inversion procedure.

Here we go directly to the calibration results of the ice-thickness inversion. OGGM currently has only one free parameter to tune for the ice thickness inversion: Glen’s creep parameter A. This is very similar to [Farinotti_etal_2009], with the difference that we are not calibrating A for each glacier individually: we tried to do that, but didn’t manage (yet).

Land-terminating glaciers¶

Currently, A is varied until the glacier-wide thickness computed by OGGM minimizes the RMSD with GlaThiDa. We removed all ice-caps and marine-terminating glaciers from the dataset in order to avoid these specific cases (135 glaciers left). Here are the results of the calibration (left: volume-area scaling, right: OGGM):

We can see that OGGM has a slightly lower score than volume-area scaling (VAS). This is due to the presence of a couple of outliers. In particular, the thickest glacier in GlaThiDa and VAS is strinkingly thin in OGGM: this is the Black Rapids glacier in Canada, which is quite well mentioned in the literature because it is a surging glacier. Closer inspection in OGGM reveals that this glacier has one of the lowest mass-balance gradient. CRU precipitation is 849 mm yr \(^{1}\) (after application of the 2.5 correction factor!), which I assume is too low.

We recall here that one central difference between our approach and that of [HussFarinotti_2012] is that we use real climate data to compute the apparent mass-balance, and thus have glacier specific mass-balance gradients. This is a strength but can also become a burden. The mass-balance gradient depends mostly on precipitation, but also on temperature and its seasonal cycle. Here I show the apparent mass-balance gradient in the ablation area of all GlaThiDa glaciers:

Is the MB gradient related to the error OGGM makes in comparison to GlaThiDa?

No not really. And this is similar with all other glacier characteristics I could look at until now. The error that OGGM makes is not easily attributable to specific causes... It it would, that would be great! Indeed, this would allow maybe to define a rule for the calibration factor A. If you have ideas at which parameter to look at, let me know!

VAS vs OGGM¶

I don’t want the calibration to be too altered by the Black Rapids outlier, so I removed it from the calibration set. The plot now looks better, and after this little hack OGGM is even slightly better than VAS:

What is really interesting however is that OGGM and VAS are incredibly similar in their dissimilarity with GlaThiDa. So similar that if we plot the two approaches together on a scatter, one could argue that the thousands of lines of code of OGGM really aren’t worth the effort ;). I think that I have to talk to David Bahr about this, he will surely be pleased.

A problem with large glaciers?¶

These results for the reference glaciers availble in GlaThiDa were somehow OK, but the issue with the Black Rapids glacier made me wonder a little bit and I decided to have a closer look. In the figure below I compare the OGGM inversion results with VAS for all glaciers I have at hand for the experiment (ITMIX + WGMS + GlaThiDa, land-terminating, no ice-caps).

After many hours (days?) of searching for a reasonable explanation for this behavior, I had to renounce. I will have to make a test under controlled conditions to see if this is actually a bug in OGGM or if it is something inherent to the methodology.

Putting all this together¶

We find an optimised factor for A of 3.22 (i.e. our A is three times larger than the standard of 2.4e-24 [Cuffey_Paterson_2010]). This makes sense since we do not consider the sliding velocity in the inversion, which means that we need an ice which is less stiff to compensate.

The inversion procedure in OGGM is not designed to provide a distributed estimate of glacier thickness: in particular, the glacier widths in OGGM are not always geometrical as in Farinotti et al. (2009): they are also corrected so that the altitude-area distribution of the glacier is preserved.

We show a few examples of the normal inversion procedure in OGGM, which is meant to provide the intput to a flowline model:

Distributed ice thickness¶

In a final step, we had to somehow interpolate the flowline thicknesses to the ITMIX map. Here again we made the choice to keep as much logic as possible within the standard OGGM framework and defined a new task, oggm.tasks.distribute_thickness(). We compute the ice thickness on the OGGM local map and simply interpolate our results on the high resolution ITMIX grid only at the very end.

The distribution works as follows:

• for each pixel, the closest flowlines thicknesses within a 100m altitude range are interpolated using an inverse distance weight
• this value is then corrected with a factor depending on the distance to the glacier outlines
• finally, the thickness is corrected with a factor \(1 / \alpha^{\frac{N}{N+2}}\) as in Farinotti et al. (2009). Since I had no time to check if this correction works better, I submitted two versions of the interpolation (with or without local slope correction).

An example of the interpolation for Columbia glacier with (left) or without (right) slope correction:

For the ice caps, the interpolation methods lead to similar results. However, it is very likely that the flowline methodology used by OGGM is not working properly.

Conclusions¶

We were able to provide an estimate of ice thickness for all glaciers except Starbuck in Antarctica. I am less confident about the distributed maps than the total volume estimates. I am also much less confident about the ice-caps (they are probably totally crap) and also less confident about the larger glaciers than the smaller ones.

References¶

Entity tasks¶

Global tasks¶

Classes¶

Mass-balance¶

Mass-balance models in the oggm.core.models.massbalance module