Coding a simple temperature index model

Coding a simple temperature index model#

Today we are going to code a simple temperature index model. We will start with creating the most basic degree day model (DDM). After that we will guide you step-by-step to add more complexity to the DDM.

Also if you’re new to python and/or a beginner in programming in general, there is no need to worry. We will take this step by step. We will start simple and it will get more challenging as we go. The goal of the activity is to get a grasp of temperature index models and get some hands-on experience (not to finish the notebook: so please go at your own pace).

Like always, we first load the packages we need. (Though this can be done later, doing it now keeps the document nice and tidy.)

# imports and defaults
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
sns.set_style('ticks')
sns.set_context('notebook')

In addition to loading the packages, we will need climate data for our degree day model and mass balance data to calibrate our model. We will use the AWS data at Zhadang:

Click the link to download: aws_data_zhadang_localtime.nc. You should put the file in the same folder as the notebook! We will demonstrate how to do this.

These data are from a weather station on a glacier in Tibet (Zhadang) that Fabien maintained during his PhD.

Lets have a look what is in the file. As you can see (after executing the next cell) there are quite some variables in this file, but we will only use what we need.

Note: we gave you this data because you learnt this format and xarray in the previous notebook. If you prefer a csv file, you can download it here.

The data is has a hourly resolution. As the name of the model states, we need daily data in a degree day model. Let’s compute that:

There are more variable here than we actually need. Let’s have a look at the SR50 variable first. The sonic ranger (SR50) is a stationary device that measure the changes in surface elevation by measuring the travel return time of a pulse:

ds.SR50.plot();

../_images/15d8ce331a1359e78c249d4dcfd60b3f02b4cb2369daa4d674bedbcfdd0cddac.png

As we can see, there is quite some data gaps. This is because the SR50 was buried in snow during the first winter. Therefore, lets focus on the melting season 2012, form May 15th until August 15th:

ds = ds.sel(time=slice('2012-05-15', '2012-08-15'))

ds.TEMP.plot();
plt.title('Daily temperature'); plt.ylabel('°C');

../_images/486f25b7c034ebf7ad88f01b97412972aa893faa16c4fb80db51112cfa566c09.png

ds.SR50.plot();
plt.title('Surface height (0 = ice)'); plt.ylabel('[m]');

../_images/12a6cc997c6e59f9f001de60c4b2c2e531889e4136efa1c0cfdd18d0c903a4f9.png

This plot shows surface elevation, but what we want to simulate is melt in units of water. The usual unit of melt is mm water equivalent (or, better: kg m-2). Let us convert surface elevation into melt, using the following density values:

snow: 200 kg m-4
ice/firn: 850 kg m-4

Let us do this for you (no need to understand what is going on):

# No need to understand the code below for now
is_snow = ds.SR50.where(ds.SR50 > 0, -1) > 0
snow_melt = - ds.SR50.where(ds.SR50 > 0, 0) * 200
ice_melt = - ds.SR50.where(ds.SR50 <= 0, 0) * 850
observed_melt = snow_melt + ice_melt - snow_melt.isel(time=0)

observed_melt.plot();
plt.title('Cumulative melt'); plt.ylabel('[kg m-2]');

../_images/d8f44c55e3dee83d385237e260ddc2e0a53a2bffe5cd72c0fb1ef39bdf433c93.png

As we can see, the melting at the beginning of the season (in presence of snow) is actually slower than the melting of ice (when converted to water equivalent).

1. Simplest DDM#

For the simplest DDM, we are going to use a single factor for the entire melting period. The DDM formulation is very simple:

\[Melt = f \cdot PDD\]

Where PDD is the sum total of daily average temperatures above 0°C in a given time period and f is a melting factor. Melt is computed in mm water equivalent (w.e. or kg m-2) per day (d-1), and therefore the melt factor has units kg m-2 K-1 d-1.

Let’s start by computing the “positive degree days”, which is the sume of temperature above the melt threshold (here 0°):

# The function "where" works as follows: where our statement is true, keep the data. When false, replace by 0
pdd = ds.TEMP.where(ds.TEMP > 0, 0)

pdd.plot();
plt.title('Daily degree days'); plt.ylabel('°C');

../_images/cc6d323a7945111b07393392a236eb121494c843075bc48a3937a7d850888c81.png

As documented since quite some time, “positive degree days” are a good indicator of ice or snow melt over a day. Let’s verify this!

Let’s assume that the melt factor is 12.

Activity: compute daily modelled melt according to the formula above. Then compute the cumulative melt over the period, and plot it on the same graph as the cumulative observed melt.

Tip: you will need the .cumsum() function

melt_f = 5
modelled_melt = (melt_f * pdd).cumsum()
observed_melt.plot(label='Cumulative melt (observations)')
modelled_melt.plot(label='Cumulative melt (model)')
plt.ylabel('[kg m-2]')
plt.legend();

../_images/1d063803f8e57640bc5b9beecef6b7fecd7743cf1c19a12207d1662530cfadc3.png

Can't go any further? Click for the answer


melt_f = 12

modelled_melt = (melt_f * pdd).cumsum()

observed_melt.plot(label='Cumulative melt (observations)')

modelled_melt.plot(label='Cumulative melt (model)')

plt.ylabel('[kg m-2]')

plt.legend();

2. Calibrated DDM#

Not too bad for a model! But I’m sure we can do a little bit better. First of all, we could compute the melt factor so that observations are matched perfectly (calibration).

Activity: compute the melt factor so that melt is matched perfectly. Replot your results.

Tip: you may need the .max() function to compute the total cumulative melt, and you may also need to compute the sum of the positive degree days over the period with .sum().

# calibrate
optimal_melt_f = observed_melt.max() / pdd.sum()

modelled_melt = (optimal_melt_f * pdd).cumsum()
observed_melt.plot(label='Cumulative melt (observations)')
modelled_melt.plot(label='Cumulative melt (better model)')

plt.ylabel('[kg m-2]')
plt.legend();

../_images/6e8bda1b3ff03cf00e2a06895c69c20df0049c85f2c559c6d0ff84a728e25b8c.png

Can't go any further? Click for the answer


optimal_melt_f = observed_melt.max() / pdd.sum()

modelled_melt = (optimal_melt_f * pdd).cumsum()

observed_melt.plot(label='Cumulative melt (observations)')

modelled_melt.plot(label='Cumulative melt (better model)')

plt.ylabel('[kg m-2]')

plt.legend();

3. Optional: different factors for ice and snow#

Even with the perfectly calibrated model, the pergormance in the first half is not very good. This is because we are not distiguishing between the snow surface (which has a higher albedo) from the ice surface. Traditionally, we use different factors for ice and snow.

Optional activity: if you are a versed python coder, you might be able to code a model which is able to do that. You will need the is_snow variable we computed earlier. You may also compute the optimal parameter yourself, but here is our suggestion for them: snow_factor=6.11, ice_factor=11.22.

snow_factor.values

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[14], line 1
----> 1 snow_factor.values

NameError: name 'snow_factor' is not defined

ice_factor.values

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[15], line 1
----> 1 ice_factor.values

NameError: name 'ice_factor' is not defined

# Compute the snow melt factor
snow_melt = observed_melt[is_snow].max()
snow_factor = snow_melt / pdd[is_snow].sum()

# Compute the ice melt facto
ice_melt = observed_melt[~is_snow].max() - observed_melt[~is_snow].min()
ice_factor = ice_melt / pdd[~is_snow].sum()

# Total model
best_model_daily = snow_factor * pdd.where(is_snow, 0) + ice_factor * pdd.where(~is_snow, 0)
best_model = best_model_daily.cumsum()

# Plot
observed_melt.plot(label='Cumulative melt (observations)')
best_model.plot(label='Cumulative melt (best model)')
plt.ylabel('[kg m-2]')
plt.legend();

../_images/afcb631f831c4ac1a3ac062c83e22b06fa551e90599f78f43b46f47f845361d5.png

Can't go any further? Click for the answer


# Compute the snow melt factor

snow_melt = observed_melt[is_snow].max()

snow_factor = snow_melt / pdd[is_snow].sum()

# Compute the ice melt facto

ice_melt = observed_melt[~is_snow].max() - observed_melt[~is_snow].min()

ice_factor = ice_melt / pdd[~is_snow].sum()

# Total model

best_model_daily = snow_factor * pdd.where(is_snow, 0) + ice_factor * pdd.where(~is_snow, 0)

best_model = best_model_daily.cumsum()

# Plot

observed_melt.plot(label='Cumulative melt (observations)')

best_model.plot(label='Cumulative melt (best model)')

plt.ylabel('[kg m-2]')

plt.legend();

Concluding remarks#

Note that a full mass balance model should take into account both accumulation (snowfall) and ablation (melt). Furthermore, distributed models like OGGM need to simulate the mass-balance at different elevations, which requires other paramters (such as lapse rate). Fundamentally, however, the equations is very simple. If you are interested, here is the code of the equation in OGGM. “Easy”, right?

Temperature index models can be extended to include more processes, such as radiation or snow aging. OGGM PhD student Lilian, for example, wrote a model which takes into account the influence of aging of the snow on the melt factor. For this, you’ll need to implement a kind of “snowpack” model, where snow height is tracked over time.

If you are finished early and want to go further, have a look at this OGGM tutorial on plotting the mass balance. In two days, we will go deeper into the implementation and calibration of the models in OGGM.

Coding a simple temperature index model

Contents

Coding a simple temperature index model#

1. Simplest DDM#

2. Calibrated DDM#

3. Optional: different factors for ice and snow#

Concluding remarks#