<jupyter_start><jupyter_text>Glacier surging experiments Goals of this notebook:- The student will be able to describe the effects of glacier surging.- The student will be able to implement glacier surging in OGGM. In this notebook we are going to explore the basic functionalities of OGGM flowline model(s) and perform a simple surging experiment. For this purpose we are going to use simple and "idealized" glacier models with simple linear mass-balance profiles.  **What is glacier surge?**A small percentage of glaciers undergo short periods of faster flow. They experience a change in morphology and surface characteristics, which sometimes leads to a marked frontal advance. The speed of the ice increases up to 10 - 1000 times of the normal velocity. Record flows are observed with velocities that exceed 100m per day. Surges happen in decadal cycles and can last for 1 up to 15 years.[(Jiskoot, 2011)](https://www.springer.com/la/book/9789048126415) We start with importing all needed modules:<jupyter_code># 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

# 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>Basics We set up a simple glacier with a linear bedrock (see [getting started with flowline models](flowline_model.ipynb)) as a setting for our experiment.<jupyter_code># This is the bed rock, linearily decreasing from 3400 m
# altitude to 1400 m, in 200 steps
nx = 200
bed_h = np.linspace(3400, 1400, nx)

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

# Let's set the model grid spacing to 100 m 
map_dx = 100

# Calculate the corresponding distance along the glacier (from the top)
distance_along_glacier = np.linspace(0, nx, nx) * 0.1  # in km

# The units of widths is in "grid points", i.e. 3 grid points = 300 m
# in our case
widths = np.zeros(nx) + 3.

# 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>Then we define the mass balance model, we set the ELA to 3000 m a.s.l. and the gradient to 4 mm/m:<jupyter_code># ELA at 3000m a.s.l., gradient 4 mm m-1
mb_model = LinearMassBalance(3000, grad=4)<jupyter_output><empty_output><jupyter_text>Here we calculate a glacier that flows only because of deformation and  without sliding:<jupyter_code># Define the year until we want to simulate the glacier.
year = 1200

# Initialize the model.
model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.)
# Simulate the glacier until
model.run_until(year)
# Store the outline of the glacier as simple_glacier_h
simple_glacier_h = model.fls[-1].surface_h
# Let's plot it!
plt.plot(distance_along_glacier, simple_glacier_h, label='Deforming glacier')
plt.plot(distance_along_glacier, bed_h, label='Bedrock', ls=':', c='k')
plt.xlabel('Distance along glacier [km]')
plt.ylabel('Altitude [m]')
plt.legend();<jupyter_output><empty_output><jupyter_text>Surging glacier experiment Since surging occurs because of increased basal sliding, a surging period can be represented in the model by increasing the sliding parameter[(Jiskoot, 2011)](https://www.springer.com/la/book/9789048126415), for example by a factor of 10. Typically surging can be assumend to happen cyclically after 100 years of normal sliding and lasts over a period of 10 years. Let's model a simple glacier that only slides:<jupyter_code># Define parameters:
# Default in OGGM: Glen's creeping parameter
glen_a = cfg.PARAMS['glen_a']

# Sliding glaciers (sliding parameter fs is nonzero)
fs = 5.7e-20

# Define the year until we want to simulate the glacier.
year = 1200<jupyter_output><empty_output><jupyter_text>Initialise the model. **Note that we pass Glens creeping parameter (`glen_a`) and the sliding parameter (`fs`) this time**<jupyter_code># Initialize the model
model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.,
                      glen_a=glen_a, fs=fs)
# Simulate the glacier until 
model.run_until(year)
# Store the final results for later use
sliding_glacier_h = model.fls[-1].surface_h<jupyter_output><empty_output><jupyter_text>We can plot both the glaciers to see if there is any difference<jupyter_code># Let's plot both models.
# The simple glacier
plt.plot(distance_along_glacier, simple_glacier_h, label='Deforming glacier')
# The sliding glacier
plt.plot(distance_along_glacier, sliding_glacier_h, label='Sliding glacier')
# The bed
plt.plot(distance_along_glacier, bed_h, label='Bedrock', ls=':', c='k')
plt.xlabel('Distance along glacier [km]')
plt.ylabel('Altitude [m]')
plt.legend();<jupyter_output><empty_output><jupyter_text>Now let us simulate a surging glacier<jupyter_code># Define parameters:
# Sliding parameter in times of slow sliding:
fs = 5.7e-20

# Sliding parameter in surging periods:
fs_surge = 5.7e-20*10

# Surging period (years):
period_s = 10

# Time span between 2 surging periods (years):
t_slow = 100

# Number of surges
no_surges = 10<jupyter_output><empty_output><jupyter_text>For the surging experiment we will have to divide the simulation into different parts. The first step is to generate an array with the years for which we want to bin the simulation. During year with no surging we go in steps of ten while during the surges we take steps of one year.<jupyter_code># 10 periods of 10 years surging every 100 years
yrs = np.arange(0, t_slow + 1, 10)
for i in range(no_surges):
    # Generate 10 one year increments from the last year.
    yrs_sliding = np.arange(t_slow+1+t_slow*i+period_s*i,
                            t_slow+period_s+1+t_slow*i+period_s*i,1)
    # Append them to the no surgning years
    yrs = np.append(yrs, yrs_sliding)

    # Generate another 100 years in 10 year increments.
    yrs_normal = np.arange(t_slow+period_s+10+t_slow*i+period_s*i,
                           2*t_slow+period_s+10+t_slow*i+period_s*i, 10)
    # Append it to the list.
    yrs = np.append(yrs, yrs_normal)<jupyter_output><empty_output><jupyter_text>If it is still unclear to you which years are "surging years", take look at the array `yrs`<jupyter_code>yrs<jupyter_output><empty_output><jupyter_text>Now it is time to set up the surging experiment. Similar to the experiments with the [mass balance gradients](mas_balance_gradient.ipynb), we cant simply change the parameters in the middle of the simulation. Instead we have to re-initialise a new model every time we want change any parameters. In this case it means running model **A** for 100 years, then running model **B** for 10 year, then back to model **A** for 100 year and so on.<jupyter_code># Surging glaciers - model initialisation
model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.,
                      glen_a=glen_a, fs=fs)
# The number of steps
nsteps = len(yrs)
# Lists for saving indermidate data
length_surge_weak = []
volume_surge_weak = []
surging_glacier_h_weak = []

# We loop for the glacier evolution.
for i, yr in enumerate(yrs):
    # Simulate the glacier until the desired year.
    model.run_until(yr)
    # Save the states
    length_surge_weak.append(model.length_m)
    volume_surge_weak.append(model.volume_km3)
    
    # If we are at the first or last year, we do noting
    if i == 0 or i == (nsteps-1):
        continue
    
    # Check if the next year is a surging year
    elif (yr-yrs[i-1]) == 10 and (yrs[i+1]-yr) == 1:
        # Save the surface height before surge.
        surging_glacier_h_weak.append(model.fls[0].surface_h)
        # If it is, initialise a "new glacier" based on the 
        # old shape, but with a new sliding parameter.
        model = FlowlineModel(model.fls[0], mb_model=mb_model,
                              y0=yr, glen_a=glen_a, fs=fs_surge)

    # Check if the next year is not a surging year
    elif (yr-yrs[i-1]) == 1 and (yrs[i+1]-yr) == 10:
        # Save the surface height before surge.
        surging_glacier_h_weak.append(model.fls[0].surface_h)
        # If it is, initialise a "new glacier" based on the 
        # old shape, but with a new sliding parameter.
        model = FlowlineModel(model.fls[0], mb_model=mb_model,
                              y0=yr, glen_a=glen_a, fs=fs)<jupyter_output><empty_output><jupyter_text>In the next figures the development of the glacier (length and volume) over the years is shown:<jupyter_code># Plot it in two subplots
from matplotlib.patches import Patch
from matplotlib.lines import Line2D

fig, axs = plt.subplots(2, 1, sharex=True, figsize=(14, 8))
# Spacing
fig.subplots_adjust(hspace=0)
# Specify the xticks.
plt.xticks(np.arange(min(yrs), max(yrs)+1, 100.0))
# Plot the glacier length
axs[0].plot(yrs, length_surge_weak, color='tab:green')
# Mark surging periods with orange patches.
for i, yr in enumerate(yrs):
    if (yr-yrs[i-1]) ==  1 and (yrs[i-1+period_s]-yrs[i+period_s-2]) == 1:
        axs[0].axvspan(yrs[i-1], yrs[i-1+period_s], color='tab:orange',
                       alpha=0.3)
# Label and grid.
plt.xlabel('Years')
axs[0].set_ylabel('Length [m]')
axs[0].grid()

# Plot the glacier volume
axs[1].plot(yrs, volume_surge_weak, color='tab:green')
# Mark surging periods with orange patches.
for i, yr in enumerate(yrs):
    if (yr-yrs[i-1]) ==  1 and (yrs[i-1+period_s]-yrs[i+period_s-2]) == 1:
        axs[1].axvspan(yrs[i-1], yrs[i-1+period_s], color='tab:orange',
                       alpha=0.3)
# Ylabel and grid.
axs[1].set_ylabel('Volume [km³]')
axs[1].grid()
# Define handles for the legend
label = [Line2D(yrs, volume_surge_weak, color='tab:green',
                label='Weakly surging glacier'),
         Patch(facecolor='tab:orange', alpha=0.3,
               edgecolor='r', label='Surging period')
        ]
# Add the legend
axs[0].legend(handles=label);<jupyter_output><empty_output><jupyter_text>Now we compare the different calculated glaciers:<jupyter_code># Plot the results
f, ax = plt.subplots(figsize=(9, 6))
plt.plot(distance_along_glacier, simple_glacier_h, label='Deforming glacier')
plt.plot(distance_along_glacier, sliding_glacier_h, label='Sliding glacier')
# Plot the last surface height of the surging glacier.
plt.plot(distance_along_glacier, surging_glacier_h_weak[-1],
         label='Weakly surging glacier')

# Add the bedrock
plt.plot(distance_along_glacier, bed_h, label='Bedrock', ls=':', c='k')
plt.xlabel('Distance along glacier [km]')
plt.ylabel('Altitude [m]')
plt.legend();<jupyter_output><empty_output><jupyter_text>We can increase the sliding parameter in order to increase the magnitude of the  surges:<jupyter_code># Define parameters:
# Sliding parameter in times of normal sliding:
fs = 5.7e-20
# Sliding parameter in surging periods:
fs_surge = 5.7e-20*50

# Initialise the glacier.
model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.,
                      glen_a=glen_a, fs=fs)
# The number of steps
nsteps = len(yrs)
# Lists for saving indermidate data
length_surge_strong = []
volume_surge_strong = []
surging_glacier_h_strong = []

# We loop for the glacier evolution.
for i, yr in enumerate(yrs):
    # Simulate the glacier until the desired year.
    model.run_until(yr)
    # Save the states
    length_surge_strong.append(model.length_m)
    volume_surge_strong.append(model.volume_km3)
    
    # If we are at the first or last year, we do noting
    if i == 0 or i == (nsteps-1):
        continue
    
    # Check if the next year is a surging year
    elif (yr-yrs[i-1]) == 10 and (yrs[i+1]-yr) == 1:
        # Save the surface height before surge.
        surging_glacier_h_strong.append(model.fls[0].surface_h)
        # If it is, initialise a "new glacier" based on the 
        # old shape, but with a new sliding parameter.
        model = FlowlineModel(model.fls[0], mb_model=mb_model,
                              y0=yr, glen_a=glen_a, fs=fs_surge)

    # Check if the next year is not a surging year
    elif (yr-yrs[i-1]) == 1 and (yrs[i+1]-yr) == 10:
        # Save the surface height before surge.
        surging_glacier_h_strong.append(model.fls[0].surface_h)
        # If it is, initialise a "new glacier" based on the 
        # old shape, but with a new sliding parameter.
        model = FlowlineModel(model.fls[0], mb_model=mb_model,
                              y0=yr, glen_a=glen_a, fs=fs)
# Plot the results
f, ax = plt.subplots(figsize=(9, 6))
# The simple and sliding glacier.
plt.plot(distance_along_glacier, simple_glacier_h, label='Deforming glacier')
plt.plot(distance_along_glacier, sliding_glacier_h, label='Sliding glacier')
# The weakly surgning glacier
plt.plot(distance_along_glacier, surging_glacier_h_weak[-1],
         label='Weak surging glacier')
# The more more strongly surging glacier.
plt.plot(distance_along_glacier, surging_glacier_h_strong[-1],
         label='Strong surging glacier')
# Add the bedrock
plt.plot(distance_along_glacier, bed_h, label='Bedrock', ls=':', c='k')
plt.xlabel('Distance along glacier [km]')
plt.ylabel('Altitude [m]')
plt.legend();<jupyter_output><empty_output><jupyter_text>Compare the different glacier types. Can you explain the differences?  In the next plot the development of the weakly surging glacier and the strongly surging glacier are shown:<jupyter_code># Plot it in two subplots
fig, axs = plt.subplots(2, 1, sharex=True, figsize=(14, 8))
# Spacing
fig.subplots_adjust(hspace=0)
# xlabel
plt.xlabel('Years')
# Glacier length
axs[0].plot(yrs, length_surge_strong, color='tab:red')
axs[0].plot(yrs, length_surge_weak, color='tab:green')
# Mark surging periods
for i, yr in enumerate(yrs):
    if (yr-yrs[i-1]) == 1 and (yrs[i-1+period_s]-yrs[i+period_s-2]) == 1:
        axs[0].axvspan(yrs[i-1], yrs[i-1+period_s], color='tab:orange',
                       alpha=0.3)
# ylabel and grid.
axs[0].set_ylabel('Length [m]')
axs[0].grid()

# Glacier volume
axs[1].plot(yrs, volume_surge_strong, color='tab:red')
axs[1].plot(yrs, volume_surge_weak, color='tab:green')
# Mark surging periods
for i, yr in enumerate(yrs):
    if (yr-yrs[i-1]) ==  1 and (yrs[i-1+period_s]-yrs[i+period_s-2]) == 1:
        axs[1].axvspan(yrs[i-1], yrs[i-1+period_s], color='tab:orange',
                       alpha=0.3)
# ylabel and grid
axs[1].set_ylabel('Volume [km³]')
axs[1].grid()
# Define handles for legend.
label = [Line2D(yrs, length_surge_strong, color='tab:red',
                label='Strongly surging glacier'), 
         Line2D(yrs, length_surge_weak, color='tab:green',
                label='Weakly surging glacier'),
         Patch(facecolor='tab:orange', alpha=0.3, edgecolor='r',
               label='Surging period')
        ]
# Add legend.
axs[1].legend(handles=label);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10477322242028124, 'fraction_numerical': 0.020928236614511125, 'mean_word_length': 4.157643854273068, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, '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': '120902', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5535, 'n_tokens_neox': 5069, 'n_words': 1569, 'path': 'oggm-edu/surging_experiment.ipynb', 'repo_name': 'LennartKoe/oggm-edu-notebooks', 'size': '15147'}