<jupyter_start><jupyter_text>Mass balance gradients (MBG) and their influence on glacier flow In the flowline model notebook we touched briefly on the mass balance gradient of a glacier. This is what we are going to take a closer look at now.If the concept of mass balance is completely new to you, have a short read about it [here](http://www.antarcticglaciers.org/glacier-processes/introduction-glacier-mass-balance/), up to the paragraph "*So what is Glacier Mass Balance?*".In this notebook we will set up a few simple runs to further explore what characteristics of glaciers change with different mass balance gradients. We will also take a look at the topics of volume and response time. Goals of this notebook:- The student will be able to explain the terms: mass balance gradient, equilibrium state and response time.- The student will be able to compute different MBGs with OGGM.- The student will be able to compute response times for different equilibirum line altitudes. First, we have to import all needed modules:<jupyter_code># Plotting libraries and plot style
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (10,8)

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

# Scientific packages
import numpy as np
import pandas as pd

# Constants
from oggm import cfg
cfg.initialize_minimal()

# 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
# Glacier shape
from oggm.core.flowline import RectangularBedFlowline<jupyter_output><empty_output><jupyter_text>Initialising the glacierWe start by setting up a simple model with a linear bedrock (see [getting started with flowline models](flowline_model.ipynb)) to generate starting point for our experiment.<jupyter_code># define horizontal resolution of the model:
# nx: number of grid points
nx = 200
# map_dx: grid point spacing in meters
map_dx = 100
# define glacier top and bottom altitudes in meters
top = 3400
bottom = 1400
# This is the bed rock, linearily decreasing points from the top altitude
# to bottom altitude, in nx steps
bed_h = np.linspace(top, bottom, nx)

# At the begining, there is no glacier so our glacier surface
# is at the bed altitude
surface_h = bed_h

# Calculate the distance along the glacier (from the top)
distance_along_glacier = np.linspace(0, nx, nx) * map_dx * 1e-3

# Define the glacier width as we did in flowline_model
# Width in meters
initial_width = 300  

# 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<jupyter_output><empty_output><jupyter_text>With this done we can initialise the glacier flowline<jupyter_code># Define our bed
initial_flowline = RectangularBedFlowline(surface_h=surface_h, bed_h=bed_h,
                                          widths=widths, map_dx=map_dx)<jupyter_output><empty_output><jupyter_text>Changing the mass balance gradient (MBG)The MBG is defined as the change of the mass balance with altitude [¹](References). It depends strongly on the climate at the glacier [²](References).Let's take a look at the effects of the MBG by simulating a few glaciers with different gradients. In the flowline notebook our glacier had a gradient of 4 mm/m so lets add a glacier with a weaker gradient and one with a stronger gradient.<jupyter_code># Define the list of MBGs we want to compare
# We will calculate models with the MBGs: 0.3, 4 and 15 mm/m. 
# These numbers can be found for real glaciers. 
gradients = [0.3, 4, 15]

# We want to run the model until year 300.
year = 300
# Lists for storing intermediate steps of the models.
mb_models = []
annual_mb_list = []
# Here the "final" models will be saved
models = []
# Colors for the graphs
colors = ['C1', 'C3', 'C5'] 
# Set the equilibrium line altitude (ELA)
ELA = 3000

for gradient in gradients:
    # Calculate the mass balance model for each gradient 
    mb_model = LinearMassBalance(ELA, grad=gradient)
    # Save it in the list
    mb_models.append(mb_model)
    # Calculate the annual mass balance along the glacier profile
    annual_mb = mb_model.get_annual_mb(surface_h) * cfg.SEC_IN_YEAR
    # Save it in the list.
    annual_mb_list.append(annual_mb)
    # The model require the initial glacier bed, a mass balance model,
    # and an initial time (the year y0)
    model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.)
    # Run the model until year.
    model.run_until(year)
    # Save the model state.
    models.append(model)
# Plot the results
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(14, 6), sharey=True)

# Plot the annual mass balance
for k, ann_mb in enumerate(annual_mb_list):
    ax1.plot(ann_mb, bed_h, color=colors[k],
             label=f'Mass balance, grad = {gradients[k]}')
# Add the ELA, where mass balance = 0
ax1.axhline(y=ELA, color='k', linestyle='--', linewidth=0.8,
            label='Equilibrium line altitude')
# Add a vertical line for the zero mass balance.
ax1.axvline(x=0, color='k', linestyle='--', linewidth=0.8)
ax1.set_xlabel('Annual mass balance (m yr-1)')
ax1.set_ylabel('Altitude (m)')
ax1.legend(loc='best');

# Plot the glaciers
# Plot the initial conditions first.
ax2.plot(distance_along_glacier, initial_flowline.surface_h,
         label='Initial glacier')
# Get the modelled flowline (model.fls[0]) and plot its new surface
for k, model in enumerate(models):
    ax2.plot(distance_along_glacier, model.fls[0].surface_h,
             label=f'Glacier after {model.yr} years. grad: {gradients[k]}',
             color=colors[k])
# Add the ELA:
ax2.axhline(y=ELA, color='k', linestyle='--', linewidth=0.8,
            label='Equilibrium line altitude')
# Add the bedrock:
ax2.plot(distance_along_glacier, bed_h, color='k', label='Bedrock',
         linestyle=':', linewidth=1.5);
ax2.set_xlabel('Distance along glacier (km)')
ax2.legend(loc='lower left');<jupyter_output><empty_output><jupyter_text>A stronger mass balance gradient implies a larger change of the mass balance with altitude. We can see this in the left plot: The annual mass balance hardly changes with altitude for the glacier with the weakest mass balance gradient (orange). On the other hand, there is a considerable difference between the top and bottom annual mass balance for the glacier with the strongest mass balance gradient (brown).This in turn affects the growth of the glacier. A strong mass balance gradient implies that more ice is added to the accumulation zone each year compared to glaciers with a weaker gradient, and the speed, and reach, of the moving ice is because of this greater. This is why the glacier with the strongest gradient exhibits the largest growth during our simulations (Brown glacier in the right plot).    What do you think: where do we find glaciers with high MBGs? Click for detailsYou will find a short explanation in this [paragraph](http://www.antarcticglaciers.org/modern-glaciers/introduction-glacier-mass-balance/SECTION_3) on AntarcticGlaciers.org.  Equilibrium state  Glaciers change their geometry as a way to adapt to the climate[³](References). If the accumulation increases, the glacier will grow further down below the ELA to increase the ablation. Similarly, if the temperature rises and the ablation increases, the glacier length will decrease.If the climate remains constant for a long enough time, glaciers will reach an equilibrium state with its climate, where the accumulation = ablation [⁴](References).With this in mind, we will take a look at how fast our glaciers, with different gradients, reach this state and compare their shapes:<jupyter_code># First we run the models until equilibrium.
for model in models:
    model.run_until_equilibrium(rate=0.006)
# Plot it
plt.figure()
# Plot the initial conditions first:
plt.plot(distance_along_glacier, initial_flowline.surface_h,
         label='Initial glacier')

for k, model in enumerate(models):
    # Get the modelled flowline (model.fls[0]) and plot its new surface
    plt.plot(distance_along_glacier, model.fls[0].surface_h, 
             label=f'Glacier with grad={gradients[k]} in eq. at year {model.yr}', 
             color=colors[k])
plt.title('Equilibrium state of glacier models')
# Add the bedrock:
plt.plot(distance_along_glacier, bed_h, ls=':', c='k',
        label='Bedrock')
plt.ylabel('Altitude (m)')
plt.xlabel('Distance along glacier (km)')
plt.legend();<jupyter_output><empty_output><jupyter_text>The different glaciers reach their equilibrium state after a different number of years. What does the figure show us? Which glacier is the thickest and longest? Let's look at specific numbers:<jupyter_code># Create lists with different properties of the glacier models
volume = []
area = []
length = []
year = []

for model in models:
    year.append(model.yr)
    volume.append(model.volume_km3)
    length.append(model.length_m)
    area.append(model.area_km2)

# Create a table with the properties of the glacier models
dic = {'gradient': gradients, 'year': year, 'volume': volume, 'area': area,
       'length': length}
table = pd.DataFrame.from_dict(dic)
table<jupyter_output><empty_output><jupyter_text>The glacier with the strongest gradient reaches the equilibrium state first. This is also the largest and longest glacier.  VolumeLet's take a look at the volume of the glaciers in relation to their gradient:<jupyter_code># Plot the volume depending on the gradient
plt.figure(figsize=[9,6])
plt.plot(table['gradient'], table['volume'])
plt.xlabel('gradient')
plt.ylabel('volume [km³]');<jupyter_output><empty_output><jupyter_text>The plot above show somewhat of a trend but it is rather coarse. Let's redo the experiment but with more glaciers. Be patient, the calculation will take a little bit longer.<jupyter_code># We will calculate models with the MBGs: 1 to 19 with stepsize 2
gradients = list(range(1, 21, 2))
# In the following lists the intermediate steps of the models will be saved
mb_models = []
annual_mb_list = []
# here the "final" models will be saved
models = []

for gradient in gradients:
    mb_model = LinearMassBalance(ELA, grad=gradient)
    mb_models.append(mb_model)
    # Calculate the annual mass balance along the glacier profile
    annual_mb = mb_model.get_annual_mb(surface_h) * cfg.SEC_IN_YEAR
    annual_mb_list.append(annual_mb)
    # The models require the initial glacier bed, a mass balance model,
    # and an initial time (the year y0)
    model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.)
    # Run the model until equilbrium
    model.run_until_equilibrium(rate=0.006)
    # Save it for later.
    models.append(model)<jupyter_output><empty_output><jupyter_text>Then create the table again<jupyter_code># create lists with different properties of the glacier models
volume = []
area = []
length = []
year = []

for model in models:
    year.append(model.yr)
    volume.append(model.volume_km3)
    length.append(model.length_m)
    area.append(model.area_km2)

dic = {'gradient': gradients, 'year': year, 'volume': volume, 'area': area, 'length': length}
table = pd.DataFrame.from_dict(dic)
table
# Plot it
plt.figure(figsize=[14,8])
plt.plot(table['gradient'], table['volume'])
plt.xlabel('MBG')
plt.ylabel('Volume [km³]');<jupyter_output><empty_output><jupyter_text>This shows us that the volume of a glacier at the equilibrium state strongly depends on the mass balance gradient.  Response timeThe glacier response time is the period of time a glacier needs to adjust its geometry to changes in mass balance caused by climate change and reach a new equilibrium state. There are a some slightly different definitions for its calculation. In this notebook we will use the definition according to Oerlemans (formula below) [⁵](References).We will now compare the response times of glaciers with different MBGs. Let's simulate a change in climate by changing the equilibrium line altitude for the three glacier simulated before:<jupyter_code># Initialization of the models used before:
gradients = [0.3, 4, 15]
# Lists for saving the glacier state
mb_models = []
#  ... and the models
models = []
# equilibrium line altitude (ELA)
ELA = 3000
for gradient in gradients:
    # Calculate the mass balance model for each gradient 
    mb_model = LinearMassBalance(ELA, grad=gradient)
    # Add it to the list.
    mb_models.append(mb_model)
    # The models require the initial glacier bed, a mass balance model,
    # and an initial time (the year y0)
    model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.)
    # Run the glacier to a first equilbrium.
    model.run_until_equilibrium(rate=0.006)
    # Save it.
    models.append(model)<jupyter_output><empty_output><jupyter_text>We now have the initial equilibrium state of our glaciers. Let's change the equilibrium line altitude in order to simulate a change in climate. In this case we move it down.            What does it mean for the climate by moving the ELA downwards? Click for a hint        By lowering the ELA we simulate a cooling of the climate, i.e. a larger part of the glacier experiences a positive annual mass balance.<jupyter_code># First we change the mass balance model, ELA decreases 
new_ELA = 2800<jupyter_output><empty_output><jupyter_text>Then we continue to simulate the glaciers with a the new ELA. To do this, we have to do a bit of a trick since we cant just change the ELA of the old model. Instead we initialise a new model for each glacier, but its flowline is based on the already simulated glaciers that reached equilibrium for the old ELA.<jupyter_code># Store the perturbed mb models
mb_models_pert = []
# Store the perturbed models
models_pert = []
# Loop over the different gradients.
for i, gradient in enumerate(gradients):
    # Calculate the mass balance model.
    mb_model = LinearMassBalance(new_ELA, grad=gradient)
    # Save it
    mb_models_pert.append(mb_model)
    # Initialise the new model with the eq. state flowline.
    model = FlowlineModel(models[i].fls[0], mb_model=mb_model,
                          y0=models[i].yr)
    # And run it until eq.
    model.run_until_equilibrium(rate=0.006)
    # Save it.
    models_pert.append(model)<jupyter_output><empty_output><jupyter_text>We can take a quick look at the glaciers at the new equilibrium state (also including the old equilibrium state)<jupyter_code># Plot it
plt.figure()
# Plot the initial conditions first:
# plt.plot(distance_along_glacier, initial_flowline.surface_h,
#          label='Initial glacier')

# Plot the old eq. state.
for k, model in enumerate(models):
    # Get the modelled flowline (model.fls[0]) and plot its surface.
    plt.plot(distance_along_glacier, model.fls[0].surface_h, 
             label=f'Initial eq. state with grad={gradients[k]} at year {model.yr}', 
             color=colors[k], ls='--')
# Plot the new eq. state.
for k, model in enumerate(models_pert):
    # Get the modelled flowline (model.fls[0]) and plot its surface.
    plt.plot(distance_along_glacier, model.fls[0].surface_h, 
             label=f'New eq. state with grad={gradients[k]} at year {model.yr}', 
             color=colors[k])

plt.title('Equilibrium state of glaciers')
# Add the bedrock:
plt.plot(distance_along_glacier, bed_h, ls=':', c='k',
        label='Bedrock')
# Add the ELA
plt.axhline(ELA, ls='--', zorder=0, label='Initial ELA')
plt.axhline(new_ELA, ls='--', zorder=0, label='New ELA', c='C2')
# Labels and legend.
plt.ylabel('Altitude (m)')
plt.xlabel('Distance along glacier (km)')
plt.legend(loc='lower left');<jupyter_output><empty_output><jupyter_text>The next step is to calculate the response times for our glaciers. One could think that it is as simple as looking at the years above in the plot a do a simple subtraction. Unfortunately this is not the case! In reality the rate at which a glacier changes is ever decreasing and  a complete equilibrium state is never really achieved. Because of this the response time is considered the time it has taken the glacier for most of the adjustment, all but a factor of $1/e$.  For numerical models like our glaciers it is common to use the [volume response time](https://www.tandfonline.com/doi/pdf/10.1080/00040851.1997.12003238?needAccess=true), from Oerlemans [⁵](References):$$\tau = t \left(V=V_2 - \frac{V_2 - V_1}{e}\right)$$where $V_1$ and $V_2$ corresponds to the glacier volume at the initial and new equilibrium state respectively. To compute the response time for our glaciers we will again have to simulate them to the new equilibrium state, but this time also save the intermediate steps. Here we do a little trick: We use the previously computed equilibrium states for the new ELA as a guess for how long we have to simulate the glaciers and then split this up into chunks and then use the `run_until(year)` function.<jupyter_code># Save the response times in a list.
response_times = []
# Loop over the different gradients again.
for i, gradient in enumerate(gradients):
    # Years, incremental steps of 5.
    yrs = np.arange(0, models_pert[i].yr, 5)
    # Number of steps we'll take.
    nsteps = len(yrs)
    # Save the intermediate steps.
    lengths = np.zeros(nsteps)
    volumes = np.zeros(nsteps)
    years = np.zeros(nsteps)
    # Initialize the glacier model again.
    mb_model = LinearMassBalance(new_ELA, grad=gradient)
    model = FlowlineModel(models[i].fls[0], mb_model, y0=models[i].yr)
    # Then we loop over the years.
    for j, year in enumerate(yrs):
        # Run the model until year.
        model.run_until(year)
        # Save the states
        years[j] = model.yr
        lengths[j] = model.length_m
        volumes[j] = model.volume_km3
        
    # The final volume difference is then
    v2 = model.volume_km3
    # This is the e-folding volume difference.
    v2_diff = v2 - (model.volume_km3 - models[i].volume_km3) / np.e
    # We want to find the year where the volume is closest to the v2_diff.
    # Take the difference
    all_vol_diff = np.abs(volumes - v2_diff)
    # Search for the smallest difference and get the index
    index = np.argmin(all_vol_diff)
    # Retreive the year
    response_time = years[index] - models[i].yr
    # Save it
    response_times.append(response_time)<jupyter_output><empty_output><jupyter_text>There is quite a lot going on in the cell above, but essentially we are simulating each glacier and saving the intermediate steps in order to find the year when most of the change in volume has occurred. Let's take a look at the response times (Remember that glaciers the are ordered from weak to strong MBG)<jupyter_code>response_times<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08248467075446547, 'fraction_numerical': 0.01018395094641429, 'mean_word_length': 4.343589743589743, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 1, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '242755', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5850, 'n_tokens_neox': 5509, 'n_words': 2245, 'path': 'oggm-edu/mass_balance_gradients.ipynb', 'repo_name': 'LennartKoe/oggm-edu-notebooks', 'size': '18755'}