<jupyter_start><jupyter_text>OGGM-Edu - Jupyter notebooks to learn and teach about glaciers   Authors**<NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>**- Author1 = {"name": "<NAME>", "affiliation": "University of Innsbruck", "email": "<EMAIL>", "orcid": "https://orcid.org/0000-0002-3211-506X"}- Author2 = {"name": "<NAME>", "affiliation": "University of Innsbruck", "email": "<EMAIL>", "orcid": ""}- Author3 = {"name": "<NAME>", "affiliation": "University of Innsbruck", "email": "<EMAIL>", "orcid": "https://orcid.org/0000-0003-0430-7950"}- Author4 = {"name": "<NAME>", "affiliation": "Georgia Tech", "email": "<EMAIL>", "orcid": "https://orcid.org/0000-0002-8780-3089"}- Author5 = {"name": "<NAME>", "affiliation": "University of Innsbruck", "email": "<EMAIL>", "orcid": ""}- Author6 = {"name": "<NAME>", "affiliation": "University of Innsbruck", "email": "<EMAIL>", "orcid": ""}- Author7 = {"name": "<NAME>", "affiliation": "University of Innsbruck", "email": "<EMAIL>", "orcid": "https://orcid.org/0000-0003-2101-2070"} Table of Contents1&nbsp;&nbsp;OGGM-Edu - Jupyter notebooks to learn and teach about glaciers1.1&nbsp;&nbsp;Authors1.2&nbsp;&nbsp;Purpose1.3&nbsp;&nbsp;Technical contributions1.4&nbsp;&nbsp;Methodology1.5&nbsp;&nbsp;Results1.6&nbsp;&nbsp;Funding1.7&nbsp;&nbsp;Keywords1.8&nbsp;&nbsp;Citation1.9&nbsp;&nbsp;Suggested next steps1.10&nbsp;&nbsp;Acknowledgements2&nbsp;&nbsp;Setup3&nbsp;&nbsp;Part 1: idealized glacier experiments3.1&nbsp;&nbsp;Glacier flowlines3.2&nbsp;&nbsp;Glacier mass balance3.3&nbsp;&nbsp;A first run3.4&nbsp;&nbsp;Experiment: influence of the mass-balance gradient on glacier size3.5&nbsp;&nbsp;Exercise: influence of glacier slope4&nbsp;&nbsp;Part 2: real-world glacier experiments4.1&nbsp;&nbsp;Glaciers as water resources: idealized climate4.1.1&nbsp;&nbsp;Setting the scene: glacier runoff and "peak water"4.1.2&nbsp;&nbsp;Setup4.1.3&nbsp;&nbsp;Define the glacier we will play with4.1.4&nbsp;&nbsp;Preparing the glacier data4.1.5&nbsp;&nbsp;Interactive glacier map4.1.6&nbsp;&nbsp;Generating a glacier in equilibrium with climate4.1.7&nbsp;&nbsp;Climate change simulation4.1.7.1&nbsp;&nbsp;Annual runoff4.1.7.2&nbsp;&nbsp;Monthly runoff4.1.8&nbsp;&nbsp;Wrapping up4.2&nbsp;&nbsp;Glaciers as water resources: projection with CMIP54.2.1&nbsp;&nbsp;Climate downscaling4.2.2&nbsp;&nbsp;Projection run4.3&nbsp;&nbsp;Part 2: take home points5&nbsp;&nbsp;References  PurposeOGGM-Edu (http://edu.oggm.org) is an educational platform providing instructors with tools that allow the demonstration of the governing principles of glacier behaviour. It offers online interactive apps and Jupyter notebooks that can be run in a browser via MyBinder or JupyterHub. This notebook demonstrates how an open source numerical glacier model ([OGGM](http://oggm.org)) can be used to tackle scientific and socially-relevant research questions in class assignments and student projects. It is built on two independent numerical experiments:1. the first section uses idealized glaciers (with controlled boundary conditions) to illustrate how climate and topography influence glacier behavior, and2. the second section relies on realistic glaciers to illustrate how downstream freshwater resources will be affected by glacier change in the future. The notebook contains code sections that serve as templates, followed by research questions and suggestions that invite the student to re-run the experiment with updated simulation parameters. This notebook (as well as all other similar notebooks featured on the OGGM-Edu platform) is meant to be adapted by teachers to match the topics of a specific lecture or workshop. Interested readers will find more example notebooks [on OGGM-Edu](https://edu.oggm.org/en/latest/notebooks_index.html).OGGM-Edu's notebooks can help students to overcome barriers posed by numerical modelling and coding, in a fun, question-driven and societally relevant way. Technical contributions- teach about fundamental concepts in glaciology in a modern, fun and interactive way- a numerical model used in research is used and abstracted (simplified) for teaching purposes- open and reproducible: the notebooks are all licensed with BSD3, and they can be run in openly available docker environments or on MyBinder MethodologyThe notebook relies on a number of standard packages of the scientific python ecosytem: [NumPy](https://numpy.org/), [SciPy](https://www.scipy.org/), [xarray](http://xarray.pydata.org/en/stable/index.html), [Pandas](https://pandas.pydata.org/), etc. The main tool used for the modelling of glaciers is the Open Global Glacier Model (OGGM, http://oggm.org, https://docs.oggm.org). We also rely on a small number of helper functions (mostly to hide some of the OGGM internals to the readers) from the OGGM-Edu python package. For more details about OGGM-Edu, visit https://edu.oggm.org. ResultsThe main outcome of this example notebook is a template, that can be used in glaciology workshops and classes. Example applications include a three-day workshop in Peru for university students without prior experience in glaciology ([repository](https://github.com/ehultee/CdeC-glaciologia), [more info](https://oggm.org/2019/12/06/OGGM-Edu-AGU/)) or the University of Michigan CLaSP 474 [Ice and Climate notebooks](https://github.com/skachuck/clasp474_w2021). Funding- Award1 = {"agency": "University of Innsbruck", "award_code": "Förderkreis 1669 – Wissen schafft Gesell schaft", "award_URL": "https://www.uibk.ac.at/foerderkreis1669"}- Award2 = {"agency": "Bundesministerium für Bildung und Forschung - BMBF", "award_code": "FKZ 01LS1602A", "award_URL": "https://www.bmbf.de"} Keywordskeywords=["glaciers", "climate", "education", "python", "numerical model"] Citation<NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>.. OGGM-Edu - Jupyter notebooks to learn and teach about glaciers. Jupyter Notebook accessed 07.05.2021 at https://github.com/OGGM/EarthCube2021 Suggested next stepsVisit https://edu.oggm.org and try the other notebooks there! Acknowledgements Thanks to the OGGM and OGGM-Edu developpers, and all the open-source tools used in this notebook. Most notably: [NumPy](https://numpy.org/), [SciPy](https://www.scipy.org/), [matplotlib](https://matplotlib.org/), [Pandas](https://pandas.pydata.org/), [GeoPandas](https://geopandas.org/), [xarray](http://xarray.pydata.org/en/stable/index.html)... Thanks to EarthCube for this opportunity.  Setup<jupyter_code># import oggm-edu helper package
import oggm_edu as edu

# Plotting libraries and plot style
import matplotlib.pyplot as plt
%matplotlib inline

import seaborn as sns
sns.set_context('notebook')
sns.set_style('ticks')

# Scientific packages
import numpy as np
import pandas as pd
import xarray as xr

# OGGM flowline model
from oggm.core.flowline import RectangularBedFlowline

# OGGM mass-balance model
from oggm.core.massbalance import LinearMassBalance

# There are several numerical implementations in OGGM core. We use the "FluxBasedModel"
from oggm.core.flowline import FluxBasedModel as FlowlineModel
# Constants
import oggm.cfg
# Initialize OGGM and do not write on screen
oggm.cfg.initialize_minimal(logging_level='CRITICAL')
# Prevent numerical inaccuracies
oggm.cfg.PARAMS['cfl_number'] = 5e-3
# Store only the relevant variables here
oggm.cfg.PARAMS['store_diagnostic_variables'] = ['volume', 'area', 'length']<jupyter_output><empty_output><jupyter_text>Part 1: Idealized glacier experiments This part of the notebook shows how idealized glacier geometries can be used to explore basic principles of ice flow and glacier behavior.This section is inspired from the OGGM-Edu notebooks "Glacier flowline modelling", "Ice flow parameters" and "Mass-balance gradient".  Glacier flowlines OGGM is a so called flowline model, which means that the glacier ice flow is assumed to happen along a representative “1.5D” flowline, as in the image below. “1.5D” here is used to emphasize that although glacier ice can flow only in one direction along the flowline, each point of the glacier has a geometrical width. This width means that flowline models are able to match the observed area-elevation distribution of true glaciers, and can parametrize the changes in glacier width with thickness changes. The first step for the experiment is now to define our idealized glacier geometry with OGGM:<jupyter_code># define horizontal resolution of the model:
# nx: number of grid points
# map_dx: grid point spacing in meters
nx = 200
map_dx = 100

# define glacier top and bottom altitudes in meters
top = 3400
bottom = 1400

# create a linear bedrock profile from top to bottom
bed_h, surface_h = edu.define_linear_bed(top, bottom, nx)

# glacier width in meters
initial_width = 600

# calculate the distance from the top to the bottom of the glacier in km
distance_along_glacier = edu.distance_along_glacier(nx, map_dx)

# Now describe the widths in units of "grid points" for the model, based on grid point spacing map_dx
widths = np.zeros(nx) + initial_width / map_dx

# Define our glacier flowine
init_flowline = RectangularBedFlowline(surface_h=surface_h, bed_h=bed_h, widths=widths, map_dx=map_dx)<jupyter_output><empty_output><jupyter_text>We created the geometry of our glacier. It is very simple: ice will flow along a constant slope, and the glacier bed shape is constant as well ("U-shaped", called [rectangular](http://docs.oggm.org/en/stable/ice-dynamics.htmlrectangular) in the model).<jupyter_code># plot the glacier bedrock profile and the initial glacier surface
plt.plot(distance_along_glacier, surface_h, label='Initial glacier')
edu.plot_xz_bed(x=distance_along_glacier, bed=bed_h);<jupyter_output><empty_output><jupyter_text>The `init_flowline` variable now contains all geometrical information needed by the model. It can give access to some attributes, which are quite useless for a non-existing glacier:<jupyter_code>print('Glacier length:', init_flowline.length_m)
print('Glacier area:', init_flowline.area_km2)
print('Glacier volume:', init_flowline.volume_km3)<jupyter_output><empty_output><jupyter_text>Glacier mass balance For the next step we will need a mass-balance model. To simplify things here, glaciers will gain mass via solid precipitation and loose mass via melt. To learn more about mass-balance processes, visit our [mass-balance notebook](https://oggm.org/oggm-edu-notebooks/oggm-edu/accumulation_and_ablation.html).To simplify things even more, we will use a simple linear function to describe the mass balance, defined by the [equilibrium line altitude](http://www.antarcticglaciers.org/glacier-processes/mass-balance/introduction-glacier-mass-balance/) (ELA) and the gradient (mass-balance change with altitude, unit [mm yr$^{-1}$ m$^{-1}$]):<jupyter_code># ELA at 3000m a.s.l., gradient 4 mm yr -1 m-1
ELA = 3000  # equilibrium line altitude in meters above sea level
altgrad = 4  # altitude gradient in mm/m
mb_model = LinearMassBalance(ELA, grad=altgrad)<jupyter_output><empty_output><jupyter_text>The OGGM mass-balance model computes the mass balance for any given altitude (unit: meters of ice per time [m s$^{-1}$], a simplification for the ice dynamics model). Lets compute the *annual* mass balance along the glacier profile:<jupyter_code>annual_mb = mb_model.get_annual_mb(surface_h) * oggm.cfg.SEC_IN_YEAR

# Plot it
plt.plot(annual_mb, bed_h, color='C2', label='Mass balance')
plt.xlabel('Annual mass balance (m yr-1)')
plt.ylabel('Altitude (m)')
plt.legend(loc='best')

# Display equilibrium line altitude, where annual mass balance = 0
plt.axvline(x=0, color='k', linestyle='--', linewidth=0.8)
plt.axhline(y=mb_model.ela_h, color='k', linestyle='--', linewidth=0.8);
plt.text(-7, ELA+30, 'ELA' )<jupyter_output><empty_output><jupyter_text>These values are not too far from observed values of accumulation and melt in many regions of the world.  A first run  To let our glacier grow in different climate and topography conditions we need a glacier evolution model. A glacier evolution model is simply a combination of the initial bed topography (the flowline) and the climate (the mass-balance model). Let's build it!<jupyter_code># The model requires the initial glacier bed, a mass-balance model, and an initial time (the year y0)
glacier_model = FlowlineModel(init_flowline, mb_model=mb_model, y0=0.)
# Let's run the model for a number of years (here 600), and store the output
_, ds = glacier_model.run_until_and_store(600)<jupyter_output><empty_output><jupyter_text>The model output is stored in an [xarray](http://xarray.pydata.org) dataset. Fortunately, it makes exploring the output relatively easy:*Tip: click on the variables and the icons on the right to explore their attributes*<jupyter_code>ds<jupyter_output><empty_output><jupyter_text>Now let's plot the glacier attributes length, area and volume over time:<jupyter_code>f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))
ds.length_m.plot(ax=ax1);
ax1.set_xlabel('Years'); ax1.set_ylabel('Length (m)');
(ds.area_m2 * 1e-6) .plot(ax=ax2);
ax2.set_xlabel('Years'); ax2.set_ylabel('Area (km$^2$)');
(ds.volume_m3 * 1e-9) .plot(ax=ax3);
ax3.set_xlabel('Years'); ax3.set_ylabel('Volume (km$^3$)');<jupyter_output><empty_output><jupyter_text>The first year of the simulation sees a big step increase in length and area: this is simply because the mass balance is positive at the top of the glacier, and we go from no glacier at all to a thin, but already 4 km long glacier. After this, the glacier constantly grows until reaching an equilibrium. This happens when the mass gain at the top of the glacier is compensated by mass loss at the tongue.Do you see the small steps in length and area (specifically after the year 400)? This is because the length of a glacier in OGGM is defined by grid points of a specific grid size, i.e the step size.What does a glacier at equilibrium look like? Let's have a look:<jupyter_code># Get the glacier geometry at the end of the simulation
geometry = glacier_model.get_diagnostics()<jupyter_output><empty_output><jupyter_text>The output of `get_diagnostics()` is a [Pandas](https://pandas.pydata.org/) data frame, which is very descriptive. We won't have a look at all output variables today, but here is an overview:<jupyter_code>geometry.head()
geometry.bed_h.plot(color='k', linestyle=':', label='Bedrock')
geometry.surface_h.plot(label='Glacier surface');
plt.legend(); plt.xlabel('Distance along flowline [m]'); plt.ylabel('Altitude [m]');<jupyter_output><empty_output><jupyter_text>What does it mean when a glacier is in equilibrium with its climate? It means that its volume doesn't change with time (which we can verify in the plots above), and that the total mass balance (averaged over the glacier area) is approximately equal to zero:<jupyter_code>mb_model.get_specific_mb(fls=glacier_model.fls)  # units of mm water equivalent per year<jupyter_output><empty_output><jupyter_text>Close enough! We are in the following situation, where accumulation at the top compensates ablation at the bottom:*Source: OGGM-Edu's [open glacier graphics](http://edu.oggm.org/en/latest/glacier_basics.html), made by <NAME> ([<NAME>ros yeux](http://atelierlesgrosyeux.com/))* Now, **a very common misconception** is that a glacier which doesn't advance or retreat also doesn't move. We are going to show that our glacier is well and alive, and does move quite a lot even if it's not growing or shrinking!<jupyter_code># plot the model glacier
plt.figure(figsize=(18, 5))
ax = plt.subplot(121)
velocity = geometry.ice_velocity * oggm.cfg.SEC_IN_YEAR
velocity.plot(ax=ax);
plt.legend(); plt.xlabel('Distance along flowline [m]'); plt.ylabel('Altitude [m]');
plt.title('Ice velocity along the glacier (m yr$^{-1}$)');
plt.subplot(122)
edu.plot_glacier_graphics(num='06', title='Ice flow compensates for mass balance')<jupyter_output><empty_output><jupyter_text>What is happening here? If enough ice piles up, it will move under the influence of gravity - very much like a fluid! At equilibrium, this ice flow compensates exactly for the mass gain (it removes ice at the top) and the mass loss (it replaces ice melt at the bottom).A few meters per year is not enough for us to notice, but satellites can:   *Source: [FYFD](https://fyfluiddynamics.com/2019/04/the-movement-of-glaciers-is-driven-by-gravity-the/)*  Experiment: Influence of the mass-balance gradient on glacier size  There are glaciers in many mountain areas of the world. To verify this by yourselves, we recommend to spend a few minutes on the [World Glacier Explorer](https://edu.oggm.org/en/latest/explorer.html), one of OGGM-Edu's web applications.In the app, we learned that glaciers can be located in maritime climates (regions with high amounts of precipitation and relatively warm) and in dry climates (rather dry and cold). These climate features can be translated to our simple framework by the mass-balance gradient:- A strong mass-balance gradient represents high amounts of solid precipitation and high amounts of melt (**maritime climate**)- A weak mass-balance gradient represents low amounts of solid precipitation and low amounts of melt (**dry climate**)Now let's translate this into our model code:<jupyter_code># Define the mass-balance gradients (MBG) we want to compare
# We will calculate models with the MBGs: 1, 4 and 15 mm yr-1 m-1.
# These values are realistic for real-world glaciers
gradients = [1, 4, 15]

# We use lists to store the output
mb_profiles = []
glacier_diagnostics = []
glacier_evolution = []

# equilibrium line altitude (ELA)
ELA = 3000

# Loop over the chosen gradients
for gradient in gradients:
    # Preparation of the mass-balance model for each gradient
    mb_model = LinearMassBalance(ELA, grad=gradient)

    # Run the simulation
    glacier_model = FlowlineModel(init_flowline, mb_model=mb_model, y0=0.)
    _, ds = glacier_model.run_until_and_store(1200)

    # Store the output for the plots later
    mb_profiles.append(mb_model.get_annual_mb(surface_h) * oggm.cfg.SEC_IN_YEAR)
    glacier_diagnostics.append(glacier_model.get_diagnostics())
    glacier_evolution.append(ds)
# Don't worry, this can take some time!<jupyter_output><empty_output><jupyter_text>This time, we ran the model for 1200 years instead of 600 - we will get back to this soon. Let's plot our glaciers first:<jupyter_code># Plot the results
f, (ax1, ax2) = plt.subplots(1, 2, figsize=[14, 6], sharey=True)

# Annual mass balance
ax = ax1
for gradient, profile in zip(gradients, mb_profiles):
    ax.plot(profile, bed_h, label='Mass balance, grad={}'.format(gradient))
# Add ELA, where mass balance = 0:
ax.axhline(y=ELA, color='k', linestyle='--', linewidth=0.8, label='Equilibrium line altitude')
ax.axvline(x=0, color='k', linestyle='--', linewidth=0.8)
ax.set_xlabel('Annual mass balance (m yr-1)'); ax.set_ylabel('Altitude (m)'); ax.legend(loc='best');

# Glacier profiles
ax = ax2
ax.plot(glacier_diagnostics[0].index, glacier_diagnostics[0].bed_h, color='k', label='Bedrock')
for gradient, diagnostics in zip(gradients, glacier_diagnostics):
    ax.plot(diagnostics.index, diagnostics.surface_h, label='Glacier surface, grad={}'.format(gradient));
ax.set_xlabel('Distance along flowline [m]'); ax.set_ylabel('Altitude [m]'); ax.legend(loc='best');
ax.axhline(y=ELA, color='k', linestyle='--', linewidth=0.8, label='Equilibrium line altitude');<jupyter_output><empty_output><jupyter_text>**Questions to address:**- Visually check the amount of melt at the terminus of the glacier. Compare this with melt values you can find in the literature.- Why are all glaciers symmetrical around the equilibrium line altitude? In the jargon, the ratio of the area above the equilibrium line altitude (ELA) divided by the entire area is called the Accumulation Area Ratio (AAR): here, it is 50%. Compare this with values you can find in the literature, and explain which of the assumptions we used can lead to this rather unusual AAR.- Now look at the top of the glacier in the accumulation area. Based on your observations, explain what the "mass-balance elevation feedback" might be, and what role it plays in our case.- Try to make a list of the mechanisms responsible for the different glacier sizes. **Bonus: Temporal evolution of the glaciers** Now let's plot the volume evolution of all three glaciers with time:<jupyter_code>for gradient, ds in zip(gradients, glacier_evolution):
    (ds.volume_m3 * 1e-9).plot(label='Glacier volume, grad={}'.format(gradient));
plt.legend(); plt.xlabel('Years'); plt.ylabel('Volume (km$^3$)');<jupyter_output><empty_output><jupyter_text>**Questions to address:**- Compare how long each glacier needs to reach equilibrium. - Why does one glacier respond much faster compared to another?- Now plot the ice velocity of all three glaciers at equilibrium. Discuss.<jupyter_code>for gradient, geometry in zip(gradients, glacier_diagnostics):
    (geometry.ice_velocity * oggm.cfg.SEC_IN_YEAR).plot(label='Velocity, grad={}'.format(gradient));
plt.legend(); plt.xlabel('Years'); plt.ylabel('Ice velocity (m yr$^{-1}$)');<jupyter_output><empty_output><jupyter_text>Exercise: Influence of glacier slope  Now repeat the analyses by varying another parameter: glacier slope.**Questions to address:**- What is the influence of the glacier slope on how large the glacier will be in equilibrium and how fast the glacier responds. Can you explain this?- Plot the ice velocities. This might help to understand what happens. *For instructors: see the code below.*<jupyter_code># Reset the mass-balance model to defaults. ELA at 3000m a.s.l., gradient 4 mm yr -1 m-1
ELA = 3000  # equilibrium line altitude in meters above sea level
altgrad = 4  # altitude gradient in mm/m
mb_model = LinearMassBalance(ELA, grad=altgrad)

# We use lists to store the output
mb_profiles = []
glacier_diagnostics = []
glacier_evolution = []

# To change the slope, we don't change anything to our bedrock parameters but the grid point resolution
factors = [0.7, 1, 1.5]
for factor in factors:
    # Define our glacier flowine
    init_flowline = RectangularBedFlowline(surface_h=surface_h, bed_h=bed_h, widths=widths, map_dx=map_dx * factor)
    # The glacier evolution model and the run
    glacier_model = FlowlineModel(init_flowline, mb_model=mb_model, y0=0.)

    # Run the simulation
    glacier_model = FlowlineModel(init_flowline, mb_model=mb_model, y0=0.)
    _, ds = glacier_model.run_until_and_store(1200)

    # Store the output for the plots later
    mb_profiles.append(mb_model.get_annual_mb(surface_h) * oggm.cfg.SEC_IN_YEAR)
    glacier_diagnostics.append(glacier_model.get_diagnostics())
    glacier_evolution.append(ds)
# Glacier profiles
f, ax = plt.subplots()
for factor, diagnostics in zip(factors, glacier_diagnostics):
    diagnostics.bed_h.plot(ax=ax, color='k', linestyle=':', label='')
    diagnostics.surface_h.plot(ax=ax, label='Glacier surface, factor={}'.format(factor));
ax.set_xlabel('Distance along flowline [m]'); ax.set_ylabel('Altitude [m]'); ax.legend(loc='best');
ax.axhline(y=ELA, color='k', linestyle='--', linewidth=0.8, label='Equilibrium line altitude');
for factor, ds in zip(factors, glacier_evolution):
    (ds.volume_m3 * 1e-9).plot(label='Glacier volume, factor={}'.format(factor));
plt.legend(); plt.xlabel('Years'); plt.ylabel('Volume (km$^3$)');
for factor, geometry in zip(factors, glacier_diagnostics):
    (geometry.ice_velocity * oggm.cfg.SEC_IN_YEAR).plot(label='Velocity, factor={}'.format(factor));
plt.legend(); plt.xlabel('Years'); plt.ylabel('Ice velocity (m yr$^{-1}$)');<jupyter_output><empty_output><jupyter_text>**More questions that could be asked:**- What do you think are other factors that could influence how fast a glacier responds?- Repeat the analysis with the glacier length and compare it to the volume. What responds faster?  Part 2: Real-world glacier experiments Part 1 explained how to use OGGM to answer theoretical questions about glaciers, using idealized experiments: the effect of a different slope, the concept of the equilibrium line altitude, the mass-balance gradient, etc. Now, how do we use OGGM to explore real glaciers? This section will give us some insights.This section is inspired from the OGGM-Edu notebooks "Glaciers as water resources" part 1 and 2.  Glaciers as water resources: idealized climate Goals of this section:- Prepare a model run for a real world glacier.- Run simulations using different idealized climate scenarios to explore the role of glaciers as water resources.- Understand the concept of "peak water".  Setting the scene: glacier runoff and "peak water" "*If glaciers melt, there won't be water in the mountains anymore*".This is a sentence that we often hear from people we meet, or sometimes even in news articles. In fact, the role of glaciers in the hydrological cycle is more complex than that. In this notebook, we will explore this question using idealized climate scenarios applied to real glaciers! Before we continue, let's have a look at the expected contribution of glaciers to local streamflow. The chart below shows an idealized scenario where first the climate is constant (t$_0$-t$_1$) and the glacier is in equilibrium with its climate. Then a warming occurs (t$_1$-t$_2$) and the glacier begins to loose mass. This graph makes a few very important points, which we will explore together in this notebook.1. When a glacier is in equilibrium, it does not contribute to the annual runoff  since the winter uptake of water, in the form of snow, is as large as the summer release through melt.2. When the climate is warming, glaciers are losing mass. This water contributes to downstream runoff, and the runoff increases.3. If climate warms even more, glaciers will continue to lose mass and become significantly smaller. When there isn't much ice left to melt each year (or when climate stabilizes), their contribution will become smaller, until it becomes zero again.4. In the new equilibrium, the annual runoff is the same as before, but the seasonal contribution changed.We will now get back to all these points together, using OGGM! *Graphic from [Huss & Hock (2018)](https://www.nature.com/articles/s41558-017-0049-x)*  Setup<jupyter_code># One interactive plot below requires Bokeh
# The rest of the notebook works without this dependency
import holoviews as hv
hv.extension('bokeh')
import geoviews as gv
import geoviews.tile_sources as gts
import salem
from oggm import utils, workflow, tasks, graphics
from oggm_edu import run_constant_climate_with_bias
# OGGM options
oggm.cfg.initialize(logging_level='WARNING')
oggm.cfg.PATHS['working_dir'] = utils.gettempdir(dirname='WaterResources')
oggm.cfg.PARAMS['min_ice_thick_for_length'] = 1  # a glacier is when ice thicker than 1m<jupyter_output><empty_output><jupyter_text>Define the glacier we will play with For this notebook we use the Hintereisferner, Austria. Some other possibilities to play with:- Hintereisferner, Austria: RGI60-11.00897 - Artesonraju, Peru: RGI60-16.02444- Mera, Nepal: RGI60-15.03591- Parlung No. 94, China: RGI60-15.11693And virtually any glacier you can find the RGI Id from, e.g. in the [GLIMS viewer](https://www.glims.org/maps/glims)! Large glaciers may need longer simulations to see changes though. For less uncertain calibration parameters, we also recommend to pick one of the many reference glaciers [in this list](https://github.com/OGGM/oggm-sample-data/blob/master/wgms/rgi_wgms_links_20200415.csv), where we make sure that observations of mass balance are better matched. Let's start with Hintereisferner first and you'll be invited to try out your favorite glacier at the end of this notebook.<jupyter_code># Hintereisferner
rgi_id = 'RGI60-11.00897'<jupyter_output><empty_output><jupyter_text>Preparing the glacier data This can take up to a few minutes on the first call because of the download of the required data:<jupyter_code># We pick the elevation-bands glaciers because they run a bit faster - but they create more step changes in the area outputs
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/elev_bands/qc3/pcp2.5/no_match'
gdir = workflow.init_glacier_directories([rgi_id], from_prepro_level=5, prepro_border=80, prepro_base_url=base_url)[0]<jupyter_output><empty_output><jupyter_text>Interactive glacier map  A first glimpse on the glacier of interest.*Tip: You can use the mouse to pan and zoom in the map*<jupyter_code>sh = salem.transform_geopandas(gdir.read_shapefile('outlines'))
(gv.Polygons(sh).opts(fill_color=None, color_index=None) * gts.tile_sources['EsriImagery'] *
 gts.tile_sources['StamenLabels']).opts(width=800, height=500, active_tools=['pan', 'wheel_zoom'])<jupyter_output><empty_output><jupyter_text>For OGGM, glaciers are "1.5" dimensional along their flowline:<jupyter_code>fls = gdir.read_pickle('model_flowlines')
graphics.plot_modeloutput_section(fls);<jupyter_output><empty_output><jupyter_text>Generating a glacier in equilibrium with climate  Let's prepare a run with the `run_constant_climate_with_bias` tasks from the `oggm_edu` package. It allows us to run idealized temperature and precipitation correction scenarios in an easy way. First, let's decide on a temperature evolution:<jupyter_code>years = np.arange(400)
temp_bias_ts = pd.Series(years * 0. - 2, index=years)
temp_bias_ts.plot(); plt.xlabel('Year'); plt.ylabel('Temperature bias (°C)');<jupyter_output><empty_output><jupyter_text>Not much to see here! The `temp_bias_ts` variable describes a temperature bias that will be applied to the standard climate (see below). Here the bias is -2° all along because we want to run a so-called "spinup" run, to let the glacier grow and make sure that our glacier is in dynamical equilibrium with its climate at the end of the simulation. Let's go:<jupyter_code># file identifier where the model output is saved
file_id = '_spinup'

# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs
tasks.run_with_hydro(gdir,  # Run on the selected glacier
                     temp_bias_ts=temp_bias_ts,  # the temperature bias to apply to the average climate
                     run_task=run_constant_climate_with_bias,  # which climate scenario? See following notebook for other examples
                     y0=2009, halfsize=10,  # Period which we will average and constantly repeat
                     store_monthly_hydro=True,  # Monthly ouptuts provide additional information
                     output_filesuffix=file_id);  # an identifier for the output file, to read it later<jupyter_output><empty_output><jupyter_text>OK, so there is quite some new material in the cell above. Let's focus on the most important points:- We run the model for 400 years (as defined by our control temperature time series)- The model runs with a constant climate averaged over 21 years (2 times `halfsize` + 1) for the period 1999-2019- We apply a cold bias of -2°C. Indeed, Hintereisferner is in strong disequilibrium with the current climate and such a bias is needed to make the glacier grow (more on this later) Let's have a look at the output now:<jupyter_code>with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix=file_id)) as ds:
    # The last step of hydrological output is NaN (we can't compute it for this year)
    ds = ds.isel(time=slice(0, -1)).load()<jupyter_output><empty_output><jupyter_text>There are plenty of variables in this dataset! We can list them with:<jupyter_code>ds<jupyter_output><empty_output><jupyter_text>*Tip: you can click on a variable and show it's attribute with the "page" button on the right.*The `time` and `month_2d` variables are coordinates, and the other variables are either provided as additional information (e.g. `calendar_month`, we will get back to this), or they are providing the actual data. For instance, we can plot the annual evolution of the volume and length of our glacier:<jupyter_code>fig, axs = plt.subplots(nrows=2, figsize=(10, 7), sharex=True)
ds.volume_m3.plot(ax=axs[0]);
ds.length_m.plot(ax=axs[1]);
axs[0].set_xlabel(''); axs[0].set_title(f'{rgi_id}'); axs[1].set_xlabel('Years');<jupyter_output><empty_output><jupyter_text>The glacier grows and almost doubles its volume. After 400 years, it is in equilibrium.This spinup simulation won't be analyzed further: let's jump to the heart of the topic: **what happens with glacier runoff when climate is warming**?  Climate change simulation We start by generating a temperature change scenario. We define a short period of constant temperature (at -2°C, to be consistent with above), followed by a linear increase (2.5 degrees in about 150 years) and a stabilization:<jupyter_code>ny_s = 50  # Start
ny_t = 150  # Trend
ny_e = 102  # Stabilisation
temp_bias_ts = np.concatenate([np.full(ny_s, -2.), np.linspace(-2, 0.5, ny_t), np.full(ny_e, 0.5)])
temp_bias_ts = pd.Series(temp_bias_ts, index=np.arange(ny_s + ny_t + ny_e))
temp_bias_ts.plot(); plt.xlabel('Year'); plt.ylabel('Temperature bias (°C)');
# file identifier where the model output is saved
file_id = '_lin_temp'

# We are using the task run_with_hydro to store hydrological outputs along with the usual glaciological outputs
tasks.run_with_hydro(gdir,
                     temp_bias_ts=temp_bias_ts,  # the temperature bias timeseries we just created
                     run_task=run_constant_climate_with_bias,  # which climate scenario? See following notebook for other examples
                     y0=2009, halfsize=10,  # Period which we will average and constantly repeat
                     store_monthly_hydro=True,  # Monthly ouptuts provide additional information
                     init_model_filesuffix='_spinup',  # We want to start from the glacier in equibrium we created earlier
                     output_filesuffix=file_id);  # an identifier for the output file, to read it later<jupyter_output><empty_output><jupyter_text>Let's have a look at our glacier evolution:<jupyter_code>with xr.open_dataset(gdir.get_filepath('model_diagnostics', filesuffix=file_id)) as ds:
    # The last step of hydrological output is NaN (we can't compute it for this year)
    ds = ds.isel(time=slice(0, -1)).load()
fig, axs = plt.subplots(nrows=3, figsize=(10, 7), sharex=True)
ds.volume_m3.plot(ax=axs[0]);
ds.length_m.plot(ax=axs[1]);
temp_bias_ts.plot(ax=axs[2], c='C3');
axs[0].set_xlabel(''); axs[0].set_title(f'{rgi_id}'); axs[1].set_xlabel(''); axs[2].set_xlabel('Years'); axs[2].set_ylabel('Temperature bias (°C)');<jupyter_output><empty_output><jupyter_text>In this scenario, the glacier would melt almost entirely (it looses about 90% of its volume). Note also that the glacier continues to adjust after the temperature has stabilized. **What are the implications for downstream runoff?**  Annual runoff As glaciers retreat, they contribute to sea-level rise (visit the [World Glaciers Explorer](https://edu.oggm.org/en/latest/explorer.html) OGGM-Edu app for more information!). This is however not what we are interested in for this notebook. Instead we will explore an important local impact of glaciers: their impact on streamflow.  Let's take a look at some of the hydrological outputs computed by OGGM. We start by creating a pandas DataFrame of all "1D" (annual) variables in the output dataset:<jupyter_code>sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]
df_annual = ds[sel_vars].to_dataframe()<jupyter_output><empty_output><jupyter_text>Then we can select the hydrological variables and sum them to get the total annual runoff:*Note that we convert the annual runoff from kilos to mega-tonnes (Mt)*<jupyter_code># Select only the runoff variables
runoff_vars = ['melt_off_glacier', 'melt_on_glacier', 'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
# Convert them to megatonnes (instead of kg)
df_runoff = df_annual[runoff_vars] * 1e-9
# We smooth the output, which is otherwize noisy because of area discretization
df_runoff = df_runoff.rolling(31, center=True, min_periods=1).mean()
fig, ax = plt.subplots(figsize=(10, 3.5), sharex=True)
df_runoff.sum(axis=1).plot(ax=ax);
plt.ylabel('Runoff (Mt)'); plt.xlabel('Years'); plt.title(f'Total annual runoff for {rgi_id}');<jupyter_output><empty_output><jupyter_text>The hydrological variables are computed on the largest possible area that was covered by glacier ice during the simulation. This is equivalent to the runoff that would be measured at a fixed-gauge hydrological station at the glacier terminus.The total annual runoff consists of the following components:- Melt off-glacier: snow melt on areas that are now glacier free (i.e. 0 in the year of largest glacier extent, in this example at the start of the simulation)- Melt on-glacier: ice + seasonal snow melt on the glacier- Liquid precipitation on- and off-glacier (the latter being zero at the year of largest glacial extent, in this example at the start of the simulation)<jupyter_code>f, ax = plt.subplots(figsize=(10, 6));
df_runoff.plot.area(ax=ax, color=sns.color_palette('rocket')); plt.xlabel('Years');
plt.ylabel('Runoff (Mt)'); plt.title(rgi_id);<jupyter_output><empty_output><jupyter_text>The glacier length and volume decrease from year 50 onwards until about year 250 - this is the **glacier retreat** phase. Afterwards, length and volume stabilize at a constant value indicating that the glacier has **reached equilibrium**. Now study the graph above and it's four main components.**Questions to address:**- When is "peak water" reached? Does this have anything to do with a change in the temperature trend itself? If not, what is the reason for "peak water" to occur?- Verify that the total fixed-gauge runoff (remember what this means?) is the same at the beginning and the end of the simulation. Why is that so? What changed in between?- The temperature stabilizes at year 200, yet all hydrological variables continue to change after that. What is happening there?- What is the contribution of liquid precipitation at the beginning and the end of the simulation? What changed?- Can you verify that the net glacier contribution to runoff (i.e. - $\Delta M / \Delta t$ is zero at the beginning and the end of the simulation and positive in between? Why do we expect this behavior?*Hint: compare the results to the [idealised runoff graphic that we introduced at the beginning of this chapter](runoff-graphic)*<jupyter_code># Solution of the last question above
glacier_mass = ds.volume_m3.to_series() * oggm.cfg.PARAMS['ice_density'] * 1e-9  # In Megatonnes, Mt
glacier_mass = glacier_mass.rolling(31, center=True, min_periods=1).mean()
(- glacier_mass.diff()).plot()
plt.axhline(y=0, color='k', ls=':')
plt.ylabel('Annual glacier mass change (Mt yr$^{-1}$)')
plt.xlabel('Years'); plt.title('Glacier contribution to annual runoff');<jupyter_output><empty_output><jupyter_text>Monthly runoff The "2D" variables contain the same hydrological data, but at monthly resolution (dimensions [time, month]). For example, monthly runoff can be computed as:<jupyter_code># Select only the runoff variables and convert them to megatonnes (instead of kg)
monthly_runoff = (ds['melt_off_glacier_monthly'] + ds['melt_on_glacier_monthly'] +
                  ds['liq_prcp_off_glacier_monthly'] + ds['liq_prcp_on_glacier_monthly'])
monthly_runoff = monthly_runoff.rolling(time=31, center=True, min_periods=1).mean() * 1e-9
monthly_runoff.clip(0).plot(cmap='Blues', cbar_kwargs={'label': 'Runoff (Mt)'}); plt.xlabel('Months'); plt.ylabel('Years');<jupyter_output><empty_output><jupyter_text>But be aware, something is a bit wrong with this plot (some of you may have noticed: maximum melt happens in the fall, which is unusual). That's because the OGGM coordinates are hydrological months - a side effect of using a so called mass balance year, which begins in October and ends in September (i.e. encompasses the entire accumulation and ablation season). We can correct for this by shifting the months forward:<jupyter_code># This should work in both hemispheres maybe?
ds_roll = ds.roll(month_2d=ds['calendar_month_2d'].data[0] - 1, roll_coords=True)
ds_roll['month_2d'] = ds_roll['calendar_month_2d']

# Select only the runoff variables and convert them to megatonnes (instead of kg)
monthly_runoff = (ds_roll['melt_off_glacier_monthly'] + ds_roll['melt_on_glacier_monthly'] +
                  ds_roll['liq_prcp_off_glacier_monthly'] + ds_roll['liq_prcp_on_glacier_monthly'])
monthly_runoff = monthly_runoff.rolling(time=31, center=True, min_periods=1).mean() * 1e-9
monthly_runoff.clip(0).plot(cmap='Blues', cbar_kwargs={'label': 'Runoff (Mt)'}); plt.xlabel('Months'); plt.ylabel('Years');<jupyter_output><empty_output><jupyter_text>As we can see, the runoff is approximately zero during the winter months, while relatively high during the summer months. Now let's compare the actual runoff to total precipitation over the basin:<jupyter_code># Compute total precipitation (Snow + Liquid)
tot_precip = (ds_roll['liq_prcp_off_glacier_monthly'] + ds_roll['liq_prcp_on_glacier_monthly'] +
              ds_roll['snowfall_off_glacier_monthly'] + ds_roll['snowfall_on_glacier_monthly'])
tot_precip *= 1e-9  # in Mt
# Plot these data at year 0
yr = 0
r = monthly_runoff.sel(time=yr)
p = tot_precip.sel(time=yr)

f, ax = plt.subplots(figsize=(10, 6));
r.plot(ax=ax, color='C3', label='Monthly runoff', linewidth=3);
p.plot(ax=ax, color='C0', label='Monthly precipitation', linewidth=3);
ax.fill_between(r.month_2d, r, p, where=(p >= r), facecolor='C0', interpolate=True, alpha=0.5)
ax.fill_between(r.month_2d, r, p, where=(r > p), facecolor='C3', interpolate=True, alpha=0.5)
plt.ylabel('Mt yr$^{-1}$'); plt.legend(loc='best');
plt.xlabel('Month'); plt.title(f'Total monthly runoff and precipitation at year {yr}');<jupyter_output><empty_output><jupyter_text>At the beginning of the simulation, the glacier is in equilibrium and its annual contribution to runoff is zero. In fact, this means that the blue area in the plot above is equal to the red area: **the glaciated basin releases water during the summer, that accumulated in form of snow in the winter months**.In this region (Ötztal Alps, Austria), precipitation is relatively constant throughout the year. **Discuss what the implications could be in other climates of the world, in particular in climates with a strong seasonality of precipitation.** Now let's have a look at the seasonal change in runoff with time:<jupyter_code>f, ax = plt.subplots(figsize=(10, 6));
cmap = sns.color_palette('magma', 3)
for i, yr in enumerate([0, 120, 300]):
    monthly_runoff.sel(time=yr).plot(ax=ax, color=cmap[i], label=f'Year {yr}')
plt.ylabel('Mt yr$^{-1}$'); plt.legend(loc='best');
plt.xlabel('Month'); plt.title('Total monthly runoff change with time');<jupyter_output><empty_output><jupyter_text>**Questions to explore:**- Verify that the annual runoff is the same at year 0 than at year 300. What changed in between?- What is responsible for the switch from summer melt to spring melt?- How does the runoff change at "peak water" (approx. year 120)?- Discuss the implications of this shift in seasonality for Hintereisferner and for other regions of the world.**Now repeat this notebook (ideally by copying it first) with other regions of the world, where climate might be very different.**Note that our explanations are very glacier-specific. So if you change the glacier by using another *rgi_id*, you also might need to adapt the explanations!**You can also change the parameters of the simulations!**What happens when the temperature change is larger/faster, smaller/slower? What are the implications for "peak water"?  Wrapping up Now let's conclude this notebook by reproducing now this famous plot by [<NAME> (2018)](https://www.nature.com/articles/s41558-017-0049-x) with our data.*Graphic from [<NAME> (2018)](https://www.nature.com/articles/s41558-017-0049-x)*<jupyter_code># Make a dataframe out of the xarray dataset for simplicity
df_monthly_runoff = pd.DataFrame(monthly_runoff.data, index=ds.time, columns=monthly_runoff.month_2d)

# Create new columns for each season (we define seasons from the perspective of the northern Hemisphere)
df_monthly_runoff['Spring (Apr-Jun)'] = df_monthly_runoff[[4, 5, 6]].sum(axis=1)
df_monthly_runoff['Summer (Jul-Sep)'] = df_monthly_runoff[[7, 8, 9]].sum(axis=1)
df_monthly_runoff['Winter (Oct-Mar)'] = df_monthly_runoff[[10, 11, 12, 1, 2, 3]].sum(axis=1)
df_monthly_runoff.columns.name = 'Season'

# Plot it
# Attention the seasons and labels might need to be adapted if a glacier in a different region is used!
f, ax = plt.subplots(figsize=(10, 6));
df_monthly_runoff[['Spring (Apr-Jun)', 'Summer (Jul-Sep)', 'Winter (Oct-Mar)']].plot.area(ax=ax, color=sns.color_palette('rocket'));
plt.xlabel('Years'); plt.ylabel('Runoff (Mt)'); plt.title('Runoff by season');<jupyter_output><empty_output><jupyter_text>And now the actual plot:<jupyter_code>f, (ax1, ax2) = plt.subplots(2, 1, gridspec_kw={'height_ratios': [2, 1]}, figsize=(10, 8), sharex=True)

p1 = df_monthly_runoff[['Spring (Apr-Jun)', 'Summer (Jul-Sep)', 'Winter (Oct-Mar)']].sum(axis=1)
p2 = df_monthly_runoff['Summer (Jul-Sep)'] # July to September is here the melt season
p1.plot(ax=ax1, color='C0', label='Annual runoff', linewidth=3)
p2.plot(ax=ax1, color='C3', label='Melt-season runoff', linewidth=3);
ax1.hlines([p1.loc[0], p2.loc[0]], 0, 300, color=['C0', 'C3'], linestyle=[':']);
ax1.legend(); ax1.set_xlabel(''); ax1.set_ylabel('Runoff (Mt)');

(temp_bias_ts + 2).plot(ax=ax2, color='C2', label='Temperature', linewidth=3);
ax2.legend(loc='lower right'); ax2.set_xlabel('Years'); ax2.set_ylabel('$\Delta T (°C)$');
sns.despine();<jupyter_output><empty_output><jupyter_text>Note that the two plots above are very glacier-specific. If you use a glacier with another climate, you might need to adapt the plots, because we are always interested in those months where the runoff from the glacier contributes most to the total runoff. Example given, for *Artesonraju in Peru*, we are most interested in June, July, August and September as this is the dry season, where it almost does not rain at all.   Glaciers as water resources: projection with CMIP5 You have now learned how to simulate and analyse a specific glacier under an idealized climate. We will now take this a step further and simulate two different glaciers, located in different climatic regions, forced with CMIP5 (Coupled Model Intercomparison Project) climate projections.We begin by initializing the glacier directories:<jupyter_code># We keep Hintereisferner from earlier, but also add a new glacier
rgi_ids = [rgi_id, 'RGI60-15.02420']
gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=5, prepro_border=80, prepro_base_url=base_url)<jupyter_output><empty_output><jupyter_text>`gdirs` now contain two glaciers, one in Central Europe and one in the Eastern Himlayas:<jupyter_code>gdirs<jupyter_output><empty_output><jupyter_text>We can take a quick look at the new glacier:<jupyter_code>sh = salem.transform_geopandas(gdirs[1].read_shapefile('outlines'))
(gv.Polygons(sh).opts(fill_color=None, color_index=None) * gts.tile_sources['EsriImagery'] *
 gts.tile_sources['StamenLabels']).opts(width=800, height=500, active_tools=['pan', 'wheel_zoom'])<jupyter_output><empty_output><jupyter_text>Climate downscaling Before we run our simulation we have to process the climate data for the glacier  i.e. downscale it: (This can take some time)<jupyter_code>from oggm.shop import gcm_climate
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))
    workflow.execute_entity_task(gcm_climate.process_cmip_data, gdirs,
                                 # Name file to recognize it later
                                 filesuffix='_CCSM4_{}'.format(rcp),
                                 # temperature projections
                                 fpath_temp=ft,
                                 # precip projections
                                 fpath_precip=fp,
                                 );<jupyter_output><empty_output><jupyter_text>Projection run With the downscaling complete, we can run the forced simulations:<jupyter_code>for rcp in ['rcp26', 'rcp45', 'rcp60', 'rcp85']:
    rid = '_CCSM4_{}'.format(rcp)
    workflow.execute_entity_task(tasks.run_with_hydro, gdirs,
                                 run_task=tasks.run_from_climate_data, ys=2020,
                                 # use gcm_data, not climate_historical
                                 climate_filename='gcm_data',
                                 # use the chosen scenario
                                 climate_input_filesuffix=rid,
                                 # this is important! Start from 2020 glacier
                                 init_model_filesuffix='_historical',
                                 # recognize the run for later
                                 output_filesuffix=rid,
                                 # add monthly diagnostics
                                 store_monthly_hydro=True);
# Create the figure
f, ax = plt.subplots(figsize=(18, 7), sharex=True)
# Loop over all scenarios
for i, rcp in enumerate(['rcp26', 'rcp45', 'rcp60', 'rcp85']):
    file_id = f'_CCSM4_{rcp}'
    # Open the data, gdirs[0] correspond to Hintereisferner.
    with xr.open_dataset(gdirs[0].get_filepath('model_diagnostics', filesuffix=file_id)) as ds:
        # Load the data into a dataframe
        ds = ds.isel(time=slice(0, -1)).load()

    # Select annual variables
    sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]
    # And create a dataframe
    df_annual = ds[sel_vars].to_dataframe()

    # Select the variables relevant for runoff.
    runoff_vars = ['melt_off_glacier', 'melt_on_glacier', 'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
    # Convert to mega tonnes instead of kg.
    df_runoff = df_annual[runoff_vars].clip(0) * 1e-9
    # Sum the variables each year "axis=1", take the 11 year rolling mean
    # and plot it.
    df_runoff.sum(axis=1).rolling(window=11).mean().plot(ax=ax, label=rcp, color=sns.color_palette("rocket")[i])

ax.set_ylabel('Annual runoff (Mt)'); ax.set_xlabel('Year'); plt.title(rgi_id); plt.legend();<jupyter_output><empty_output><jupyter_text>For Hintereisferner, runoff continues to decrease throughout the 21st-century for all scenarios, indicating that "peak water" has already been reached sometime in the past. This is the case for many European glaciers. What about our unnamed glacier in the Himalayas?<jupyter_code># Create the figure
f, ax = plt.subplots(figsize=(18, 7), sharex=True)
# Loop over all scenarios
for i, rcp in enumerate(['rcp26', 'rcp45', 'rcp60', 'rcp85']):
    file_id = f'_CCSM4_{rcp}'
    # Open the data, gdirs[1] correspond to the unnamed glacier.
    with xr.open_dataset(gdirs[1].get_filepath('model_diagnostics', filesuffix=file_id)) as ds:
        # Load the data into a dataframe
        ds = ds.isel(time=slice(0, -1)).load()

    # Select annual variables
    sel_vars = [v for v in ds.variables if 'month_2d' not in ds[v].dims]
    # And create a dataframe
    df_annual = ds[sel_vars].to_dataframe()

    # Select the variables relevant for runoff.
    runoff_vars = ['melt_off_glacier', 'melt_on_glacier',
                   'liq_prcp_off_glacier', 'liq_prcp_on_glacier']
    # Convert to mega tonnes instead of kg.
    df_runoff = df_annual[runoff_vars].clip(0) * 1e-9
    # Sum the variables each year "axis=1", take the 11 year rolling mean
    # and plot it.
    df_runoff.sum(axis=1).rolling(window=11).mean().plot(ax=ax, label=rcp, color=sns.color_palette("rocket")[i])

ax.set_ylabel('Annual runoff (Mt)'); ax.set_xlabel('Year'); plt.title(gdirs[1].rgi_id); plt.legend();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10293466756470701, 'fraction_numerical': 0.01720388723179063, 'mean_word_length': 5.196017646357458, 'pattern_counts': {'":': 34, '<': 217, '<?xml version=': 0, '>': 219, 'https://': 35, 'lorem ipsum': 0, 'www.': 9, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '46', 'dir': 'jupyter-structured-clean-dedup', 'id': '144089', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 17025, 'n_tokens_neox': 16003, 'n_words': 5669, 'path': 'FM_01_OGGM-Edu_learn_and_teach_about_glaciers.ipynb', 'repo_name': 'pat-schmitt/EarthCube2021', 'size': '52304'}