Run OGGM with GCM data

In this example, we illustrate how to do a typical “projection run”, i.e. using GCM data (here CMIP5 and CMIP6). There are two important steps:

  • simulate all glaciers from their inventory date up to a “present day” geometry (New in version 1.4: this is already available in the pre-processed directories)

  • simulate their future evolution from the present day state

# Libs
import xarray as xr
import matplotlib.pyplot as plt

# Locals
import oggm.cfg as cfg
from oggm import utils, workflow, tasks
from oggm.core import gcm_climate
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In [1], line 8
      6 import oggm.cfg as cfg
      7 from oggm import utils, workflow, tasks
----> 8 from oggm.core import gcm_climate

ImportError: cannot import name 'gcm_climate' from 'oggm.core' (/usr/local/pyenv/versions/3.10.7/lib/python3.10/site-packages/oggm/core/__init__.py)

Pre-processed directories

Let’s do a run for two Himalayan glaciers: Ngojumba and Khumbu.

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

# Here we override some of the default parameters
# How many grid points around the glacier?
# Make it large if you expect your glaciers to grow large:
# here, 40 is more than enough
cfg.PARAMS['border'] = 40

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

# RGI glaciers: Ngojumba and Khumbu
rgi_ids = utils.get_rgi_glacier_entities(['RGI60-15.03473', 'RGI60-15.03733'])

# Go - get the pre-processed glacier directories
gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=5)

New in version 1.4: the _historical runs

As of v1.4, the level 5 files now come with a pre-computed model run from the RGI outline date to the last possible date given by the historical climate data. These files are stored in the directory with a _historical suffix. Let’s compile them for our two glaciers:

ds = utils.compile_run_output(gdirs, input_filesuffix='_historical')
ds.volume.plot(hue='rgi_id');

Each RGI glacier has an “inventory date”, the time at which the ouline is valid. It can change between glaciers:

gdirs[0].rgi_date, gdirs[1].rgi_date
ds.time[-1]

The command that was run on preprocessing was the run_from_climate_data which will run OGGM from this inventory date up to a desired date (here, the last year with valid historical climate data). So the data above is strictly equivalent to:

workflow.execute_entity_task(tasks.run_from_climate_data, gdirs, 
                             output_filesuffix='_my_spinup',  # to use the files as input later on
                            );
ds = utils.compile_run_output(gdirs, input_filesuffix='_my_spinup')
ds.volume.plot(hue='rgi_id');

One thing to remember here is that although we try hard to avoid spin-up issues, the glacier after the inversion is not in a perfect dynamical state. Some variable (in particular glacier length) might need some years to adjust:

ds.length.plot(hue='rgi_id');

Download and process GCM data from CMIP5

A typical use case for OGGM will be to use climate model output (here CMIP5). We use some of the files we store online here, but you can use whichever you want. The above mentioned storage contains information about the data, how to cite them and tabular summaries of the available GCMs.

We download the files and bias correct them to the calibration data (here: CRU, see process_gcm_data which is called by process_cmip_data). This step is very important, as the model is very sensitive to temperature variability.

bp = 'https://cluster.klima.uni-bremen.de/~oggm/cmip5-ng/pr/pr_mon_CCSM4_{}_r1i1p1_g025.nc'
bt = 'https://cluster.klima.uni-bremen.de/~oggm/cmip5-ng/tas/tas_mon_CCSM4_{}_r1i1p1_g025.nc'
for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    # Download the files
    ft = utils.file_downloader(bt.format(rcp))
    fp = utils.file_downloader(bp.format(rcp))
    # bias correct them
    workflow.execute_entity_task(gcm_climate.process_cmip_data, gdirs, 
                                 filesuffix='_CMIP5_CCSM4_{}'.format(rcp),  # recognize the climate file for later
                                 fpath_temp=ft,  # temperature projections
                                 fpath_precip=fp,  # precip projections
                                 );

Projection runs

We now run OGGM under various scenarios and start from the end year of the historical run:

for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    rid = '_CMIP5_CCSM4_{}'.format(rcp)
    workflow.execute_entity_task(tasks.run_from_climate_data, gdirs, ys=2020, 
                                 climate_filename='gcm_data',  # use gcm_data, not climate_historical
                                 climate_input_filesuffix=rid,  # use the chosen scenario
                                 init_model_filesuffix='_historical',  # this is important! Start from 2020 glacier
                                 output_filesuffix=rid,  # recognize the run for later
                                );

Plot model output

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))
for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    rid = '_CMIP5_CCSM4_{}'.format(rcp)
    ds = utils.compile_run_output(gdirs, input_filesuffix=rid)
    ds.isel(rgi_id=0).volume.plot(ax=ax1, label=rcp);
    ds.isel(rgi_id=1).volume.plot(ax=ax2, label=rcp);
plt.legend();

New: and now from CMIP6!

Very similar steps, although some links have changed and RCPs are now SSPs! See this table for a summary of available GCMs.

bp = 'https://cluster.klima.uni-bremen.de/~oggm/cmip6/GCM/CESM2/CESM2_{}_r1i1p1f1_pr.nc'
bt = 'https://cluster.klima.uni-bremen.de/~oggm/cmip6/GCM/CESM2/CESM2_{}_r1i1p1f1_tas.nc'
for ssp in ['ssp126', 'ssp585']:  # Removed 'ssp245', 'ssp370' because the files are large!
    # Download the files
    ft = utils.file_downloader(bt.format(ssp))
    fp = utils.file_downloader(bp.format(ssp))
    # bias correct them
    workflow.execute_entity_task(gcm_climate.process_cmip_data, gdirs, 
                                 filesuffix='_CMIP6_CESM2_{}'.format(ssp),  # recognize the climate file for later
                                 fpath_temp=ft,  # temperature projections
                                 fpath_precip=fp,  # precip projections
                                 );
for ssp in ['ssp126', 'ssp585']:
    rid = '_CMIP6_CESM2_{}'.format(ssp)
    workflow.execute_entity_task(tasks.run_from_climate_data, gdirs, ys=2020, 
                                 climate_filename='gcm_data',  # use gcm_data, not climate_historical
                                 climate_input_filesuffix=rid,  # use the chosen scenario
                                 init_model_filesuffix='_historical',  # this is important! Start from 2020 glacier
                                 output_filesuffix=rid,  # recognize the run for later
                                );
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))
for ssp in ['ssp126', 'ssp585']:
    rid = '_CMIP6_CESM2_{}'.format(ssp)
    ds = utils.compile_run_output(gdirs, input_filesuffix=rid)
    ds.isel(rgi_id=0).volume.plot(ax=ax1, label=ssp);
    ds.isel(rgi_id=1).volume.plot(ax=ax2, label=ssp);
plt.legend();

What’s next?