10 minutes to… a glacier change projection with GCM data#

In this example, we illustrate how to do a typical “projection run”, i.e. using GCM data. Here we will first use already bias-corrected CMIP6 data from ISIMIP3b and than show how alternatives like the original CMIP5 and CMIP6 data can be used.

There are three important steps:

  • download the OGGM pre-processed directories containing a pre-calibrated and spun-up glacier model

  • download the climate projections (and bias correct them in the case of CMIP5 or CMIP6)

  • simulate the future glacier evolution from the present day state to the end of the century (2020-2100)

Tags: beginner, projections, CMIP, workflow

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

# Locals
import oggm.cfg as cfg
from oggm import utils, workflow, tasks, DEFAULT_BASE_URL
from oggm.shop import gcm_climate

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')  # print less log messages than the default 

# 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 = ['RGI60-15.03473', 'RGI60-15.03733']

# Go - get the pre-processed glacier directories
# You have to explicitly indicate the url from where you want to start from
gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=5, prepro_base_url=DEFAULT_BASE_URL)
2024-08-25 21:07:15: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-08-25 21:07:15: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-08-25 21:07:15: oggm.cfg: Multiprocessing: using all available processors (N=4)
2024-08-25 21:07:17: oggm.workflow: init_glacier_directories from prepro level 5 on 2 glaciers.
2024-08-25 21:07:17: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 2 glaciers

The _spinup_historical runs#

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. In case of the new default climate dataset GSWP3_W5E5, this is until the end of 2019, so the volume is computed until January 1st, 2020. These files are stored in the directory with a _spinup_historical suffix (see the “10 minutes to… a dynamical spinup” tutorial for context).

Let’s compile them into a single file for our two glaciers:

ds = utils.compile_run_output(gdirs, input_filesuffix='_spinup_historical')
vol_ref2000 = ds.volume / ds.volume.sel(time=2000) * 100
vol_ref2000.plot(hue='rgi_id');
plt.ylabel('Volume (%, reference 2000)');
2024-08-25 21:07:18: oggm.utils: Applying global task compile_run_output on 2 glaciers
2024-08-25 21:07:18: oggm.utils: Applying compile_run_output on 2 gdirs.
../../_images/9f63ea03969344de315286efad8c71a3f49d8047bc880d2222ecf4faf18e76f2.png

Each RGI glacier has an “inventory date”, the time at which the ouline is valid:

gdirs[0].rgi_date, gdirs[1].rgi_date
(2000, 2002)

The glacier volume and area estimates before that date are highly uncertain and serve the purpose of spinup only! In the “10 minutes to… a dynamical spinup” tutorial, we talk about why. For now, these files are perfect for our purpose, since we plan to start our simulation in 2020.

Download and process GCM data from ISIMIP3b (bias-corrected CMIP6)#

A typical use case for OGGM will be to use climate model output (here bias-corrected CMIP6 GCMs from ISIMIP3b). We use the files we mirrored in Bremen here, but you can use whichever you want. From ISIMIP3b, we have 5 GCMs and 3 SSPs on the cluster. You can find more information on the ISIMIP website. Let’s download the data:

# you can choose one of these 5 different GCMs:
# 'gfdl-esm4_r1i1p1f1', 'mpi-esm1-2-hr_r1i1p1f1', 'mri-esm2-0_r1i1p1f1' ("low sensitivity" models, within typical ranges from AR6)
# 'ipsl-cm6a-lr_r1i1p1f1', 'ukesm1-0-ll_r1i1p1f2' ("hotter" models, especially ukesm1-0-ll)
member = 'mri-esm2-0_r1i1p1f1' 

for ssp in ['ssp126', 'ssp370','ssp585']:
    # bias correct them
    workflow.execute_entity_task(gcm_climate.process_monthly_isimip_data, gdirs, 
                                 ssp = ssp,
                                 # gcm member -> you can choose another one
                                 member=member,
                                 # recognize the climate file for later
                                 output_filesuffix=f'_ISIMIP3b_{member}_{ssp}'
                                 );
2024-08-25 21:07:18: oggm.workflow: Execute entity tasks [process_monthly_isimip_data] on 2 glaciers
2024-08-25 21:07:19: oggm.workflow: Execute entity tasks [process_monthly_isimip_data] on 2 glaciers
2024-08-25 21:07:21: oggm.workflow: Execute entity tasks [process_monthly_isimip_data] on 2 glaciers

The advantage of using ISIMIP3b data is that they have been bias-corrected by the ISIMIP consortium. Since we are using the W5E5 dataset as the baseline climate in OGGM v1.6, there is no need for us to bias correct any further. If you want to bias-correct the projections yourself or want to have a larger variety of GCMs, you can also use the original CMIP5 or CMIP6 GCMs.

If you want to know which historical data you are using, you can ask OGGM:

gdirs[0].get_climate_info()
{'baseline_climate_source': 'GSWP3_W5E5',
 'baseline_yr_0': 1901,
 'baseline_yr_1': 2019,
 'baseline_climate_ref_hgt': 5376.0,
 'baseline_climate_ref_pix_lon': 86.75,
 'baseline_climate_ref_pix_lat': 28.25}

Projection runs#

We now run OGGM under various scenarios starting from the end year of the historical spin-up run:

for ssp in ['ssp126', 'ssp370', 'ssp585']:
    rid = f'_ISIMIP3b_{member}_{ssp}'
    workflow.execute_entity_task(tasks.run_from_climate_data, gdirs,
                                 climate_filename='gcm_data',  # use gcm_data, not climate_historical
                                 climate_input_filesuffix=rid,  # use the chosen scenario
                                 init_model_filesuffix='_spinup_historical',  # this is important! Start from 2020 glacier
                                 output_filesuffix=rid,  # recognize the run for later
                                );
2024-08-25 21:07:22: oggm.workflow: Execute entity tasks [run_from_climate_data] on 2 glaciers
2024-08-25 21:07:23: oggm.workflow: Execute entity tasks [run_from_climate_data] on 2 glaciers
2024-08-25 21:07:24: oggm.workflow: Execute entity tasks [run_from_climate_data] on 2 glaciers

Plot model output#

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))
# Pick some colors for the lines
color_dict={'ssp126':'blue', 'ssp370':'orange', 'ssp585':'red'}
for ssp in ['ssp126','ssp370', 'ssp585']:
    rid = f'_ISIMIP3b_{member}_{ssp}'
    # Compile the output into one file
    ds = utils.compile_run_output(gdirs, input_filesuffix=rid)
    # Plot it
    ds.isel(rgi_id=0).volume.plot(ax=ax1, label=ssp, c=color_dict[ssp]);
    ds.isel(rgi_id=1).volume.plot(ax=ax2, label=ssp, c=color_dict[ssp]);
plt.legend();
2024-08-25 21:07:25: oggm.utils: Applying global task compile_run_output on 2 glaciers
2024-08-25 21:07:25: oggm.utils: Applying compile_run_output on 2 gdirs.
2024-08-25 21:07:25: oggm.utils: Applying global task compile_run_output on 2 glaciers
2024-08-25 21:07:25: oggm.utils: Applying compile_run_output on 2 gdirs.
2024-08-25 21:07:25: oggm.utils: Applying global task compile_run_output on 2 glaciers
2024-08-25 21:07:25: oggm.utils: Applying compile_run_output on 2 gdirs.
../../_images/4b88a2cadf245eedc4b9020fa04b82ab1dfd5e9051bd1577686ed779619b9a5b.png

That’s it! Your OGGM projection in 10 minutes.

Have 5 minutes more? Download and process GCM data from CMIP5 or CMIP6#

ISIMIP data is very useful because it is bias corrected. Furthermore, it offers daily data (which we will soon make available in OGGM).

But you may want a higher diversity of models or scenarios: for this, you may also use the CMIP5 or CMIP6 GCMs directly. These need to be bias-corrected first to the applied baseline climate (see process_gcm_data). This relatively simple bias-correction is automatically done by process_cmip_data and is very important, as the model is very sensitive to temperature variability (see the following blogpost for more details).

  • CMIP5 has 4 different RCP scenarios and a variety of GCMs, online you can find them here. The above mentioned storage contains information about the data, how to cite them and tabular summaries of the available GCMs.

  • CMIP6 has 4 different SSP scenarios, see this table for a summary of available GCMs. There are even some CMIP6 runs that go until 2300.

Note, that the CMIP5 and CMIP6 files are much larger than the ISIMIP3b files. This is because we use a simple processing trick for the ISIMIP3b GCM files as we only save the glacier gridpoints, instead of the entire globe for CMIP5 and CMIP6.0

Therefore: run the following code only if it is ok to download a few gigabytes of data. Set the variable below to true to run it.

download_cmip5_data = False  # set to True to run the code below
if download_cmip5_data:

    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'

    color_dict_rcp={'rcp26':'blue', 'rcp45':'violet', 'rcp85':'red'}

    # Download and bias correct the data
    for rcp in ['rcp26', 'rcp45', 'rcp85']: # 'rcp60' would also be available
        # 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
                                     );

    # Run OGGM
    for rcp in ['rcp26', 'rcp45',  'rcp85']: #'rcp60',
        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
    f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 4))
    for rcp in ['rcp26', 'rcp45', 'rcp85']: #'rcp60',
        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, c=color_dict_rcp[rcp]);
        ds.isel(rgi_id=1).volume.plot(ax=ax2, label=rcp, c=color_dict_rcp[rcp]);
    plt.legend();

Now, the same for CMIP6 but instead of RCPs, now SSPs and again with another GCM:

(Attention! This may take some time …) Set the variable below to true to run it.

download_cmip6_data = False  # set to True to run the code below
if download_cmip6_data:
    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'

    # Download and bias correct the data
    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, 
                                     #year_range=('1979', '2014'),
                                     filesuffix='_CMIP6_CESM2_{}'.format(ssp),  # recognize the climate file for later
                                     fpath_temp=ft,  # temperature projections
                                     fpath_precip=fp,  # precip projections
                                     );

    # Run OGGM
    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
                                    );

    # Plot
    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, c=color_dict[ssp]);
        ds.isel(rgi_id=1).volume.plot(ax=ax2, label=ssp, c=color_dict[ssp]);

    plt.legend();

What’s next?#