# Ice Dynamics, Part 2: Mass Balance and Dynamic Equilbirium

Introducing a mass balance, there can now be an equilbirium profile of the glacier, whose properties we will now explore. We’ll start by making the slope a bit larger, so that things fit a bit better on our domain.

This notebook has been adapted from https://edu.oggm.org/en/latest/notebooks_ice_flow_parameters.html and https://edu.oggm.org/en/latest/notebooks_flowline_intro.html

You can find and run the notebook from https://www.github.com/skachuck/clasp474_w2021.

# import oggm-edu helper package
import oggm_edu as edu

# import modules, constants and set plotting defaults
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (9, 6)  # Default plot size

# Scientific packages
import numpy as np

# Constants
from oggm import cfg
cfg.initialize_minimal()

# OGGM models
from oggm.core.massbalance import LinearMassBalance
from oggm.core.flowline import FluxBasedModel, RectangularBedFlowline, TrapezoidalBedFlowline, ParabolicBedFlowline

# There are several numerical implementations in OGGM core. We use the "FluxBasedModel"
FlowlineModel = FluxBasedModel

2021-03-23 18:02:06: oggm.cfg: Reading default parameters from the OGGM params.cfg configuration file.
2021-03-23 18:02:06: oggm.cfg: Multiprocessing switched OFF according to the ENV variable OGGM_USE_MULTIPROCESSING
2021-03-23 18:02:06: oggm.cfg: Multiprocessing: using slurm allocated processors (N=1)


## Basics

First we set-up a simple run with a standard sloping bedrock:

### Glacier bed

# define horizontal resolution of the model:
# nx: number of grid points
# map_dx: grid point spacing in meters
nx = 400
map_dx = 100

distance_along_glacier = edu.distance_along_glacier(nx, map_dx)

# define glacier top and bottom altitudes in meters
top = 3400
bottom = 1400
# create a linear bedrock profile from top to bottom
bed_h_slope, surface_h_slope = edu.define_linear_bed(top, bottom, nx)
initial_width = 300 #width in meters
# Now describe the widths in "grid points" for the model, based on grid point spacing map_dx
widths = np.zeros(nx) + initial_width/map_dx


### The initial mass balance

# ELA at 3000m a.s.l., gradient 4 mm m-1
ELA = 3000 #equilibrium line altitude in meters above sea level


## Running the Model

# Reinitialize the model
mb_flowline = RectangularBedFlowline(surface_h=bed_h_slope, bed_h=bed_h_slope, widths=widths, map_dx=map_dx)
model = FlowlineModel(mb_flowline, mb_model=mb_model_new, y0=0.)
# Year 0 to 1000 in 1 years step
yrs = np.arange(0, 1001, 5)
# Array to fill with data
nsteps = len(yrs)
length_mb = np.zeros(nsteps)
vol_mb = np.zeros(nsteps)
# Loop
for i, yr in enumerate(yrs):
model.run_until(yr)
length_mb[i] = model.length_m
vol_mb[i] = model.volume_km3
# I store the final results for later use
mb_glacier_h = model.fls[-1].surface_h

edu.glacier_plot(x=distance_along_glacier, bed=bed_h_slope, model=model, mb_model=mb_model_new, init_flowline=mb_flowline)

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(yrs, length_mb, label='Mass Balance');
ax1.set_xlabel('Years')
ax1.set_ylabel('Length (m)');
ax1.legend()
ax2.plot(yrs, vol_mb);
ax2.set_xlabel('Years')
ax2.set_ylabel('Volume (km3)');


At the very beginning, there’s a step change in the glacier’s length, as the glacier immediately extends to the ELA (since OGGM doesn’t distinguish between snow and ice). A period of relative stability where the thickness behind the ELA is builing up so that enough ice can advect beyond the ELA without immediately melting. Then, after about 1000 years, the length settles into it’s equilibrium state.

What happens if the altitidue gradient increases so that there’s more snowfall at higher altitudes, but more melt at lower ones? We can start a new model from where the previous one left off, but change the mass balance gradient.

Q1 What do you think will happen? Will the glacier get longer or shorter?

# ELA at 3000m a.s.l., gradient 5 mm m-1
ELA = 3000 #equilibrium line altitude in meters above sea level
# Reinitialize the model
mbchange_flowline = RectangularBedFlowline(surface_h=mb_glacier_h, bed_h=bed_h_slope, widths=widths, map_dx=map_dx)
model = FlowlineModel(mbchange_flowline, mb_model=mb_model_change, y0=0.)


Note that where mbchange_flowline is created, we’ve switched the initial surface (surface_h) to be the surface at the end of the last simulation (mb_glacier_h).

We’ll run it for a shorter time since we’re starting a bit closer to equilibrium.

# Year 0 to 500 in 5 years step
yrs = np.arange(0, 501, 5)
# Array to fill with data
nsteps = len(yrs)
length_mbchange = np.zeros(nsteps)
vol_mbchange = np.zeros(nsteps)
# Loop
for i, yr in enumerate(yrs):
model.run_until(yr)
length_mbchange[i] = model.length_m
vol_mbchange[i] = model.volume_km3
# I store the final results for later use
mbchange_glacier_h = model.fls[-1].surface_h

# Plot the final result:
plt.plot(distance_along_glacier, mb_glacier_h, label='4 mm/m')
plt.plot(distance_along_glacier, mbchange_glacier_h, label='5 mm/m', ls='--')
plt.plot(distance_along_glacier, bed_h_slope, 'k:', label='Bedrock')

plt.legend(); plt.xlabel('Distance along glacier (km)'); plt.ylabel('Altitude (m)');
plt.title('Glacier Thicknesses at Equilibrium');


Q2 Evaluate your hypothesis from Question 9.

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(yrs, length_mbchange);
ax1.set_xlabel('Years')
ax1.set_ylabel('Length (m)');
ax1.legend()
ax2.plot(yrs, vol_mbchange);
ax2.set_xlabel('Years')
ax2.set_ylabel('Volume (km3)');


Q3 Explain the shape of this plot. What happens initially? About how long does it take to re-establish equibrium (it’s response time)?

Now you change the ELA and perform the same analysis.

Note It is possible that the ELA you choose will cause the program to crash. If it does, that’s okay! Choose a different ELA, closer to the original 3000 m that you know works. If you do crash the program, try to figure out why. What happened? What does the error message say? What does it mean?

Q4 Predict what will happen to the glacier’s length and its response time

# Change the ELA!
ELA = #equilibrium line altitude in meters above sea level
# Reinitialize the model
elachange_flowline = RectangularBedFlowline(surface_h=mb_glacier_h, bed_h=bed_h_slope, widths=widths, map_dx=map_dx)
model = FlowlineModel(elachange_flowline, mb_model=mb_model_ela, y0=0.)

# Year 0 to 500 in 5 years step
yrs = np.arange(0, 501, 5)
# Array to fill with data
nsteps = len(yrs)
length_elachange = np.zeros(nsteps)
vol_elachange = np.zeros(nsteps)
# Loop
for i, yr in enumerate(yrs):
model.run_until(yr)
length_elachange[i] = model.length_m
vol_elachange[i] = model.volume_km3
# I store the final results for later use
elachange_glacier_h = model.fls[-1].surface_h

# Plot the final result:
plt.plot(distance_along_glacier, mb_glacier_h, label='ELA: 3000 m')
plt.plot(distance_along_glacier, elachange_glacier_h, label='ELA: {} m'.format(ELA), ls='--')
plt.plot(distance_along_glacier, bed_h_slope, 'k:', label='Bedrock')

plt.legend(); plt.xlabel('Distance along glacier (km)'); plt.ylabel('Altitude (m)');
plt.title('Glacier Thicknesses at Equilibrium');

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(yrs, length_elachange, label='Mass Balance');
ax1.set_xlabel('Years')
ax1.set_ylabel('Length (m)');
ax1.legend()
ax2.plot(yrs, vol_elachange);
ax2.set_xlabel('Years')
ax2.set_ylabel('Volume (km3)');


### Play

Going back to the factors of ice dynamics, change the temperature or the basal friction of the ice and reperform these experiments with the same initial mass balance and ELA, and the same perturbations. Predict, and then evaluate, what will happen to the change in glacier length and response time under the new temperature and sliding conditions.

Also, if you didn’t break the program before, try to break it now. Push its bounds and figure out what they are and why.