<jupyter_start><jupyter_text>What determines the glacial response to climate? The image below ([source](https://www.esa.int/Our_Activities/Observing_the_Earth/Glaciers_lose_nine_trillion_tonnes_of_ice_in_half_a_century)) shows changes in mountain glaciers since 1961. You can see that most glaciers lose ice (red) but there are also some that grow (green). The magnitude of the changes is highly variable!![glacier-changes](http://www.esa.int/var/esa/storage/images/esa_multimedia/images/2019/04/global_glacier_mass_loss_1961_2016/19340812-1-eng-GB/Global_glacier_mass_loss_1961_2016_node_full_image_2.jpg)The changes visible in this graph are regional and global, and are a very clear sign of climate change in the last decades. But if we look now at the change of some individual glaciers (where we have observations of the change in length), we already see that they are not the same ([source](https://glacierhub.org/2017/02/21/research-confirms-significance-climate-change-glacier-retreat/)): So why do some glaciers respond more strongly than others to the same change in climate?    There are many important factors, but one is the shape of the glacier. We have seen in the laboratory and also in [notebook 1](1_glacier_bed_slope.ipynb) glaciers of various bed slopes. In this notebook, we will see how the slope affects the response to climate change that we studied in [the previous notebook](4_a_changing_climate.ipynb). First, we import our modules as usual and configure the settings that will be the same for each glacier we investigate.<jupyter_code># The commands below are just importing the necessary modules and functions
# Plot defaults
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (9, 6)  # Default plot size
# Scientific packages
import numpy as np
import pandas as pd

# Constants
from oggm import cfg
cfg.initialize()
# OGGM models
from oggm.core.massbalance import LinearMassBalance
from oggm.core.flowline import FluxBasedModel, RectangularBedFlowline, TrapezoidalBedFlowline, ParabolicBedFlowline
# There are several solvers in OGGM core. We use the default one for this experiment
from functools import partial
FlowlineModel = partial(FluxBasedModel, min_dt=0, cfl_number=0.01)

# OGGM Edu helper functions
import oggm_edu as edu
import CdeC as cdec
# define horizontal resolution of the model:
# nx: number of grid points
# map_dx: grid point spacing in meters
nx = 200
map_dx = 100

# ELA at 3000 m a.s.l., gradient 4 mm m-1
initial_ELA = 3000 #equilibrium line altitude in meters above sea level
altgrad = 4 #altitude gradient in mm/m
mb_model = LinearMassBalance(initial_ELA, grad=altgrad)<jupyter_output><empty_output><jupyter_text>Three glaciers of differing slopes We follow similar steps as we did [before](1_glacier_bed_slope.ipynb) to configure the basic model.<jupyter_code># Define the slopes we want to compare, as in inclinacion_de_la_base
slope_tests = [0.1, 0.15, 0.2]
top = 3400 # m, the peak elevation

# Define the width of the glacier we want to study
width = 300 # m
widths = np.zeros(nx) + width / map_dx # express it for the model
# Here is where we will store the results of our test
models = []
beds = []
distances_along = []
# Colors for the graphs
colors = ['C1', 'C3', 'C5'] 

for k, slope in enumerate(slope_tests):
    # create a linear bedrock profile from top to bottom
    bottom = top - nx * map_dx * slope #m, elevation of the bottom of the incline based on the slope we defined
    bed_h, surface_h = edu.define_linear_bed(top, bottom, nx)
    beds.append(bed_h)
    # ask the model to calculate the distance from the top to the bottom of the glacier in km
    distances_along.append(edu.distance_along_glacier(nx, map_dx))
    flowline = RectangularBedFlowline(surface_h=surface_h, bed_h=bed_h, widths=widths, map_dx=map_dx)
    # The models require the initial glacier bed, a mass balance model, and an initial time (the year y0)
    models.append(FlowlineModel(flowline, mb_model=mb_model, y0=0.))
    # First, look at the state of the glacier models after the chosen amount of years
    models[k].run_until_equilibrium(rate=0.006)<jupyter_output><empty_output><jupyter_text>We plot our glaciers.<jupyter_code>f1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(21,5))
ax_list = (ax1, ax2, ax3)

for k, slope in enumerate(slope_tests):
    ax = ax_list[k]
    ax.plot(distances_along[k], models[k].fls[-1].surface_h, label='Slope = {}'.format(slope), color=colors[k])
    ax.plot(distances_along[k], beds[k], ls=':', color='k', label='Bedrock')
    ax.set_ylim(1500, 3700)<jupyter_output><empty_output><jupyter_text>Same glaciers, different climate Due to their differing slopes, our glaciers will have different shapes at equilibrium, and different times needed to reach equilibrium, for a given climate.**In the cell below, change the altitude of the ELA several times and explore how each glacier changes shape.**<jupyter_code># Setting new ELA, gradient 4 mm m-1
new_ELA = 3200 # equilibrium line altitude in meters above sea level
altgrad = 4 # altitude gradient in mm/m
new_mb_model = LinearMassBalance(new_ELA, grad=altgrad)

# Here is where we will store the results of our test
changed_models = []
lengths = []
volumes = []
# Colors for the graphs
colors = ['C1', 'C3', 'C5'] 

for k, m in enumerate(models):
    # take the glacier we already spun up
    cm = FlowlineModel(m.fls, mb_model=new_mb_model, y0=0.)
    changed_models.append(cm)
    # examine new equilibrium
    changed_models[k].run_until_equilibrium(rate=0.006) 
    lengths.append(changed_models[k].length_m)
    volumes.append(changed_models[k].volume_km3)<jupyter_output><empty_output><jupyter_text>We make a plot to show the results:<jupyter_code>f2, (ax4, ax5, ax6) = plt.subplots(1, 3, figsize=(21,5))
ax_list = (ax4, ax5, ax6)

for k, slope in enumerate(slope_tests):
    ax = ax_list[k]
    ax.plot(distances_along[k], models[k].fls[-1].surface_h, label='ELA = {}'.format(initial_ELA), color='darkorange')
    ax.axhline(y=initial_ELA, color='orange', ls=':')
    ax.plot(distances_along[k], changed_models[k].fls[-1].surface_h, label='ELA = {}'.format(new_ELA), color='darkblue')
    ax.axhline(y=new_ELA, color='blue', ls=':')
    ax.plot(distances_along[k], beds[k], ls=':', color='k', label='Bedrock')
    ax.legend(loc='best')
    ax.set_ylim(1500, 3700)
    ax.set_xlabel('Distance along glacier [km]')
    ax.set_ylabel('Elevation [m s.n.m.]')<jupyter_output><empty_output><jupyter_text>Response to a climate change **Form a hypothesis: How does a glacier's bed slope affect its response to climate change?**Which glacier will react the fastest? Which one will have the strongest response? Now we are going to expose our glaciers to sudden climate change. At first the glaciers will be in balance. After a few years (default is 50, `change_time` in cell), we will apply a sudden change in mass balance.<jupyter_code># Time
yrs = np.arange(0, 201, 5, dtype=np.float32)
nsteps = len(yrs)
change_time = 50 #when to apply the step change

# Output containers
elas = [np.zeros(nsteps) for i in range(len(slope_tests))]
lengths = [np.zeros(nsteps) for i in range(len(slope_tests))]
areas = [np.zeros(nsteps) for i in range(len(slope_tests))]
volumes = [np.zeros(nsteps) for i in range(len(slope_tests))]
water_vols = [np.zeros(nsteps) for i in range(len(slope_tests))]

# Loop
current_ELA = initial_ELA
change = -200 #m, the amount by which we want the initial ELA to change
for k, m in enumerate(models):
    ela_arr = elas[k]
    length_arr = lengths[k]
    area_arr = areas[k]
    volume_arr = volumes[k]
    water_arr = water_vols[k]
    # establish the glacier in equilibrium first
    m_perturb = FlowlineModel(m.fls, mb_model=new_mb_model, y0=0.)
    m_perturb.run_until_equilibrium(rate=0.006)
    initial_time = m_perturb.yr
    for i, yr in enumerate(yrs):
        m_perturb.run_until(initial_time + yr)
        if yr >= change_time:
            current_ELA = initial_ELA + change
        m_perturb.mb_model = LinearMassBalance(current_ELA, grad=4)
        ela_arr[i] = current_ELA
        length_arr[i] = m_perturb.length_m
        area_arr[i] = m_perturb.area_km2
        volume_arr[i] = m_perturb.volume_km3
        water_arr[i] = cdec.ice_to_freshwater(volume_arr[i])<jupyter_output><empty_output><jupyter_text>Let's graph how much water our glaciers hold during this simulation:<jupyter_code>f3, (ax7, ax8, ax9) = plt.subplots(1, 3, figsize=(21,5), sharey=True)
ax_list = (ax7, ax8, ax9)

for k, slope in enumerate(slope_tests):
    ax = ax_list[k]
    water_vol = water_vols[k]
    ax.plot(yrs, water_vol, label='Slope = {}'.format(slope), color=colors[k])
    ax.set_ylabel('Volumen de agua almacenado [L]')
    ax.set_xlabel('Año del simulación')<jupyter_output><empty_output><jupyter_text>We also make a normalized graph, to compare how much each glacier gains or loses compared to its initial volume.<jupyter_code>f4, (ax10, ax11, ax12) = plt.subplots(1, 3, figsize=(21,5), sharey=True)
ax_list = (ax10, ax11, ax12)

for k, slope in enumerate(slope_tests):
    ax = ax_list[k]
    volume = (volumes[k])/(volumes[k][0])
    ax.plot(yrs, volume, label='Slope = {}'.format(slope), color=colors[k])
    ax.set_ylabel('Size as fraction of initial volume')
    ax.set_xlabel('Year of simulation')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10273458445040215, 'fraction_numerical': 0.022412868632707774, 'mean_word_length': 4.951499680918953, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 34, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '477235', 'license': 'MIT', 'n_tokens_mistral': 3344, 'n_tokens_neox': 3128, 'n_words': 992, 'path': 'English/5_differing_reactions.ipynb', 'repo_name': 'ehultee/glacio-CdeC', 'size': '9325'}