<jupyter_start><jupyter_text>Getting started with flowline models: idealised glaciers OGGM employs what is known as a flowline model to simulate the evolution of glaciers. What this means is that the model represents the glacier as a group of lines along which the ice flows down the mountain. The flowlines connect on to, or branch off, each other forming a system along which the ice flows. In the model the flowlines are divided into finite boxes, each containing a part of the ice that makes up the glacier. The size of a box depends on how wide and deep the glacier is at the location. The model then calculates how much ice flows in to the box from boxes above and how much ice that leaves the box and goes into boxes below.The goal of this notebook is to further explain the concept of flowlines and how they are used to model glaciers. We will do this by simulating two idealised glaciers.We begin by importing the libraries we are going to need:<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')
plt.rcParams['figure.figsize'] = (7, 5)

# 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 the 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>First steps In our first experiment we want to create a glacier as simple as possible, this means a flat glacier bed with a constant slope. The glacier bedFirst we define some model parameters: the the resolution and spacing between the grid points:<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<jupyter_output><empty_output><jupyter_text>Then we define the altitude at the top and bottom of the glacier<jupyter_code># Define glacier top and bottom altitudes (in meters)
top = 3400
bottom = 1400<jupyter_output><empty_output><jupyter_text>With this information we can create the profile of our glacier bed, or the bedrock. This is simply `nx` linearly spaced altitudes between `top` and `bottom`.<jupyter_code># We want nx number of linearly spaced altitudes between the top and
# bottom altitude.
bed_h = np.linspace(top, bottom, nx)<jupyter_output><empty_output><jupyter_text>Before we have any ice and snow in our "glacier" the surface height is simply the same the height of the bed<jupyter_code>surface_h = bed_h<jupyter_output><empty_output><jupyter_text>Let's plot the bed to make sure that it looks like we expect. For this we will also calculate the distance along the glacier, a simple function of the number of grid points and the spacing between them.<jupyter_code># This gives us an array with the distances from
# the top of the glacier.
distance_along_glacier = np.linspace(0, nx, nx) * map_dx * 1e-3
# Plot the glacier bedrock profile and the initial glacier
# surface.
plt.plot(distance_along_glacier, surface_h, label='Initial glacier sfc height')
plt.plot(distance_along_glacier, bed_h, c='black', ls=':', label='Bed height')
plt.xlabel('Distance along glacier (km)')
plt.ylabel('Altitude (m)')
plt.title('Altitude to length of the glacier')
plt.legend();<jupyter_output><empty_output><jupyter_text>As we can see, the bedrock is a flat slope with a constant angle between the top and bottom of the glacier. Glacier widthThe next step is to decide how wide our glacier is and what the *shape* of the bed is. Since we're aiming for a simple glacier, we will use a rectangular "u-shaped" bed with a constant width of 300 meters. For more infor on bed shapes see the [documentation](http://docs.oggm.org/en/stable/ice-dynamics.htmlrectangular).<jupyter_code># Width in meters
initial_width = 300
# Now assign the widths to grid points which the model can use. We have
# nx grid points all of which have the value of the initial width
widths = np.zeros(nx) + initial_width
# Divide by map_dx to get the widths in "map units"
widths = widths / map_dx<jupyter_output><empty_output><jupyter_text>Glacier flowline initialisationWe can now take our bed and initialise a flowline object which OGGM uses in the simulation of the glacier.<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>The `initial_flowline` object contains all geometrical information needed by the model. We can access some of the attributes of the flowline. However since the glacier doesn't exist at the moment, they will be zero.<jupyter_code>print('Glacier length:', initial_flowline.length_m)
print('Glacier area:', initial_flowline.area_km2)
print('Glacier volume:', initial_flowline.volume_km3)<jupyter_output><empty_output><jupyter_text>Mass balanceFor the glacier to grow it needs a mass balance model. The mass balance model is responsible for adding snow and removing ice through melt on the glacier. In our case it will be a simple linear mass balance model, which is defined by the equilibrium line altitude (ELA) and an altitude gradient (in mm yr$^{-1}$ m$^{-1}$). The ELA defines at what altitude the mass balance is zero and the altitude gradient how much the mass balance changes with altitude. **More on this in upcoming notebooks!**We set the ELA of our glacier to 3000 meters and the altitude gradient to 4 mm yr$^{-1}$ m$^{-1}$.<jupyter_code># Equilibrium line altitude in meters above sea level
ELA = 3000
# Altitude gradient in mm/m
altgrad = 4
# Initialise the mass balance model.
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.  It uses the unit meters of ice per time (m s$^{-1}$), a unit which is simpler to use for the dynamical part of the model.Let us compute the *annual* mass balance along the glacier profile and plot it:<jupyter_code># Calculate the mass balance along the surface of the glacier.
# Since the get_annual_mb return the change per second we multiply with
# the number of seconds in a year.
annual_mb = mb_model.get_annual_mb(surface_h) * 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.title('Mass balance profile');<jupyter_output><empty_output><jupyter_text>Model runNow we have all the ingredients needed to run the model on our glacier. We begin by initialising the `FlowlineModel` object, which requires the initial flowline, the mass balance model and a year (Since our glacier is just starting out, we begin at year 0).<jupyter_code># The model requires 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.)<jupyter_output><empty_output><jupyter_text>Let's first run the model for one year:<jupyter_code># We want to run the model until year 1. 
year = 1
model.run_until(year)<jupyter_output><empty_output><jupyter_text>And let's plot the glacier after one year<jupyter_code># Plot the initial glacier.
plt.plot(distance_along_glacier, initial_flowline.surface_h,
         label='Initial glacier')
# Retreive the new flowline from the model and plot it. Glaciers in OGGM can
# have multiple flowlines, but currently we only have one.
current_flowline = model.fls[0]
# Plot the new surface height
plt.plot(distance_along_glacier, current_flowline.surface_h,
        label=f'Glacier after {model.yr} years')
# Plot the equilibrium line altitude
plt.axhline(mb_model.ela_h, linestyle='--', color='k', linewidth=0.8)
# Add the bedrock and axes labels:
plt.plot(distance_along_glacier, bed_h, color='k', label='Bedrock',
         linestyle=':', linewidth=1.5)
plt.xlabel('Distance along glacier (km)')
plt.ylabel('Altitude (m)')
plt.title('Altitude to length of the glacier')
plt.legend();<jupyter_output><empty_output><jupyter_text>This doesn't look much, is there even a glacier growing? We can also take a look at some of statistics of the glacier to get some more details:<jupyter_code>print('Year:', model.yr)
print('Glacier length (m):', model.length_m)
print('Glacier area (km2):', model.area_km2)
print('Glacier volume (km3):', model.volume_km3)<jupyter_output><empty_output><jupyter_text>From the statistics we can read that the glacier has a length of 4 km and covers an area of 1.2 km$^2$. The reason we can't see it in the plot above is because it is very thin. We can take a look at the thickness by first removing the bed height:<jupyter_code># Remove the bed height from the surface height and plot it along the glacier
# to show the ice thickness.
plt.plot(distance_along_glacier,
         current_flowline.surface_h - current_flowline.bed_h)
plt.xlabel('Distance along glacier (km)')
plt.ylabel('Ice thickness (m)');
plt.title('Ice thickness to distance of the glacier');<jupyter_output><empty_output><jupyter_text>Here we can see that there is thin cover of ice from the top of the glacier to 4 km down the glacier. This means that the glacier almost reaches the point where the glacier crosses the ELA (~5 km).The glacier will grow considerably in the upcoming years, and the ice thickness should become apparent even in the altitude - distance plot. Let us run the model until year 150 and take a look at the output.<jupyter_code># Run the model until year 150
year = 150
model.run_until(year)<jupyter_output><empty_output><jupyter_text>Plot the glacier:<jupyter_code># Plot the initial glacier.
plt.plot(distance_along_glacier, initial_flowline.surface_h,
         label='Initial glacier')
# Retreive the new flowline from the model and plot it. Glaciers in OGGM can
# have multiple flowlines, but currently we only have one.
current_flowline = model.fls[0]
# Plot the new surface height
plt.plot(distance_along_glacier, current_flowline.surface_h,
        label=f'Glacier after {model.yr} years')
# Plot the equilibrium line altitude
plt.axhline(mb_model.ela_h, linestyle='--', color='k', linewidth=0.8)
# Add the bedrock and axes labels:
plt.plot(distance_along_glacier, bed_h, color='k', label='Bedrock',
         linestyle=':', linewidth=1.5)
plt.xlabel('Distance along glacier (km)')
plt.ylabel('Altitude (m)')
plt.title('Altitude to length of the glacier')
plt.legend();<jupyter_output><empty_output><jupyter_text>Now we can clearly see the difference between the surface of the glacier and the bedrock/initial surface height. Let's print the same statistics about the glacier as before:<jupyter_code>print('Year:', model.yr)
print('Glacier length (m):', model.length_m)
print('Glacier area (km2):', model.area_km2)
print('Glacier volume (km3):', model.volume_km3)<jupyter_output><empty_output><jupyter_text>The glacier length and area has increased by ~20% while the volume has increased by more than 1000%. This is because the glacier has to build enough mass (or ice thickness) before it can begin to flow downhill and increase its length.Note that the model time is now at year 150. Running the model with the same input again calls the already calculated results and does not execute the method `model.run_until` again.<jupyter_code>model.run_until(150)
print('Year:', model.yr)
print('Glacier length (m):', model.length_m)<jupyter_output><empty_output><jupyter_text>If we want to simulate the glacier even longer, we have to set the specify the desired year and the model will compute the missing years.<jupyter_code># Run the model until year 500.
year = 500
model.run_until(year)
# Plot the initial glacier.
plt.plot(distance_along_glacier, initial_flowline.surface_h,
         label='Initial glacier')
# Retreive the new flowline from the model and plot it. Glaciers in OGGM can
# have multiple flowlines, but currently we only have one.
current_flowline = model.fls[0]
# Plot the new surface height
plt.plot(distance_along_glacier, current_flowline.surface_h,
        label=f'Glacier after {model.yr} years')
# Plot the equilibrium line altitude
plt.axhline(mb_model.ela_h, linestyle='--', color='k', linewidth=0.8)
# Add the bedrock and axes labels:
plt.plot(distance_along_glacier, bed_h, color='k', label='Bedrock',
         linestyle=':', linewidth=1.5)
plt.xlabel('Distance along glacier (km)')
plt.ylabel('Altitude (m)')
plt.title('Altitude to length of the glacier')
plt.legend();

print('Year:', model.yr)
print('Glacier length (m):', model.length_m)
print('Glacier area (km2):', model.area_km2)
print('Glacier volume (km3):', model.volume_km3)<jupyter_output><empty_output><jupyter_text>The glaciers has now grown considerably further down our made up mountain, well below the ELA.It is important to note that the model will not calculate back in time.Once calculated for year 500, the model will not run again for year 450 and will remain at year 500. Try running the cell below.  Does the output match what you expected?<jupyter_code>model.run_until(450)
print('Year:', model.yr)
print('Glacier length (m):', model.length_m)<jupyter_output><empty_output><jupyter_text>Starting overIt might be useful to store some intermediate steps of the evolution of the glacier for diagnostics. We do this by running the steps we did above in a loop and saving the length and volume in separate arrays. First we have to re-initialise the model:<jupyter_code># Reinitialise the model
model = FlowlineModel(initial_flowline, mb_model=mb_model, y0=0.)

# Array with year 0 to 600 in 5 year steps.
yrs = np.arange(0, 601, 5) 

# Arrays in which we store the data.
nsteps = len(yrs)
length = np.zeros(nsteps)
vol = np.zeros(nsteps)

# Loop over the years.
for i, yr in enumerate(yrs):
    model.run_until(yr)
    length[i] = model.length_m
    vol[i] = model.volume_km3

# Store the final glacier for later use
simple_glacier = model.fls[0]<jupyter_output><empty_output><jupyter_text>We can now plot the evolution of the glacier length and volume with time:<jupyter_code># Create two subplots.
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(yrs, length);
ax1.set_xlabel('Years'); ax1.set_ylabel('Length (m)');
ax2.plot(yrs, vol);
ax2.set_xlabel('Years'); ax2.set_ylabel('Volume (km3)');
f.suptitle('Time evolution of the length and volume');<jupyter_output><empty_output><jupyter_text>The glacier length exhibits a step function in the first year of simulation. This is because above the ELA, only accumulation takes places and OGGM currently does not differentiate between ice and snow. This means that all snow above the ELA from the first winter creates a bump in the length.After the first year, the glacier's length remains at first constant. This behaviour can be explained by the fact that above the ELA, the mass balance is positive, no matter how large the glacier is. In this initial stage, the ice is so thin that any flow bringing ice below the ELA will not be large enough to compensate for the high ablation rate, and any ice melts away.When the ice thickness has increased enough for the ice flow surpass the ablation rate below the ELA, the glacier length can begin to increase.<jupyter_code>'Glacier length from the top to the equilibrium line altitude ('+str(mb_model.ela_h)+' m) is: {} m'.format(length[1])<jupyter_output><empty_output><jupyter_text>After several centuries, the glacier reaches a balance with its climate. This means that its length and volume won't change anymore, as long as all physical parameters and the climate stay constant.  A first experiment  We have now seen how to setup a simple glacier and simulate it for any number of years. Now we will move a little bit closer to reality and define a glacier with changing widths. Like many real glaciers the new glacier will be wider at the top (in the accumulation area) and have a constant width below the ELA.<jupyter_code># We simply overwrite the widths we have from earlier.
# We now make our glacier 600 m wide at the top
upper_width = 600
# Convert the width in meters to width in "grid units", like before,
# and rewrite the first few points
widths[0:15] = upper_width / map_dx
# Define our new bed
wide_narrow_flowline = RectangularBedFlowline(surface_h=surface_h, bed_h=bed_h,
                                              widths=widths, map_dx=map_dx)<jupyter_output><empty_output><jupyter_text>We will now run our model with the new initial conditions (again for 600 years), and store the length and area in separate arrays for comparison:<jupyter_code># Reinitialise the model with the new input
model = FlowlineModel(wide_narrow_flowline, mb_model=mb_model, y0=0.)

# Array to fill with data
nsteps = len(yrs)
length_w = np.zeros(nsteps)
vol_w = np.zeros(nsteps)

# Loop over the years
for i, yr in enumerate(yrs):
    model.run_until(yr)
    length_w[i] = model.length_m
    vol_w[i] = model.volume_km3

# Store the final results for later use
wide_narrow_glacier = model.fls[0]<jupyter_output><empty_output><jupyter_text>Plot the simple glacier and the wide-narrow glacier to compare the results:<jupyter_code># Plot the final result:
plt.plot(distance_along_glacier, simple_glacier.surface_h,
         label='Simple glacier')
plt.plot(distance_along_glacier, wide_narrow_glacier.surface_h,
         label='Wide-narrow 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.title(f'Glacier states after {model.yr} years')
plt.legend();<jupyter_output><empty_output><jupyter_text>The glacier with a wider accumulation area is longer compared to the simple glacier area at year 600. With what you've learned so far in this notebook, can you come up with an explanation to why? Click for a hint    With a wider accumulation area the glacier mass above the ELA will increase quicker and the flow of ice to below the ELA will be larger compared to the glacier with a smaller accumulation area.    Below is the evolution of the length and volume with respect to time of the two glaciers.<jupyter_code>f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.plot(yrs, length, label='Simple glacier');
ax1.plot(yrs, length_w, label='Wide-narrow glacier');
ax1.legend(loc='best')
ax1.set_xlabel('Years')
ax1.set_ylabel('Length (m)');
ax2.plot(yrs, vol, label='Simple glacier');
ax2.plot(yrs, vol_w, label='Wide-narrow glacier');
ax2.legend(loc='best')
ax2.set_xlabel('Years')
ax2.set_ylabel('Volume (km3)');
f.suptitle('Time evolution of the length and volume of the two glaciers');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08305316242908174, 'fraction_numerical': 0.010191216642151713, 'mean_word_length': 5.349899933288859, 'pattern_counts': {'":': 0, '<': 117, '<?xml version=': 0, '>': 117, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '408419', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6127, 'n_tokens_neox': 5726, 'n_words': 2367, 'path': 'oggm-edu/flowline_model.ipynb', 'repo_name': 'LennartKoe/oggm-edu-notebooks', 'size': '19036'}