2.4. Point operations: a snow melt model

2.4.1. Precipitation and temperature

In this exercise you will use the dataset introduced in the visualisation section to build a simple spatial hydrological model, including rainfall, snowfall, snowmelt, and runoff.

Open dynamicTimeSeries.py created in the previous exercises, save it under a new name, temp.py. It reads the precipitation timeseries from disk. Add a statement that reads the temperature timeseries temp.tss from disk and assigns the temperature values (degrees celcius) to the map variable temperatureObserved. Write the results to disk, run the model, and have a look at the output.

---

Question: Which function did you use to read the temperature from disk?

1. open

2. timeinputscalar

3. timeinputnominal

4. timeinput

5. report

Correct answers: b.

Feedback:

In the dynamic (somewhere at the top), the following two lines are needed:

temperatureObserved = timeinputscalar('temp.tss',1)
self.report(temperatureObserved,'tempObs')

---

The temperature timeseries contains the temperature observed at a meteostation that is close by. However, temperature will be spatially variable, as it depends on the surface elevation. The temperature t_i at a grid cell i can be calculated as:

t_i = t_obs - l (e_i - e_obs)

In the equation, t_obs is the temperature observed at the meteostation, e_i is the elevation at a grid cell, e_obs is the elevation at the meteo station (here: 2058.1 m) and l is the temperature lapse rate (here: 0.005).

Add statements to your script temp.py to calculate a map (call it temperatureCorrection) containing the term l (e_i - e_obs). Write the map to disk and check the results.

---

Question: Where did you add the statements to calculate temperatureCorrection?

1. In the __init__.

2. In the initial.

3. In the dynamic.

4. In the initial and in the dynamic

Correct answers: b.

Feedback:

As the temperature correction remains constant over time, it needs to be calculated only once, at the start of the model run. This is why it should be put in the initial. We use a variable name starting with self, because later on the variable will be used in the dynamic. Without the self, it would be a local variable that is not available outside the initial method.

from pcraster import *
from pcraster.framework import *

class MyFirstModel(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('dem.map')

  def initial(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    elevationAboveMeteoStation = dem - elevationMeteoStation
    temperatureLapseRate = 0.005
    self.temperatureCorrection = elevationAboveMeteoStation * \
                                 temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

  def dynamic(self):
    precipitation = timeinputscalar('precip.tss',1)
    self.report(precipitation,'pFromTss')
    temperatureObserved = timeinputscalar('temp.tss',1)
    self.report(temperatureObserved,'tempObs')
    
nrOfTimeSteps=181
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

---

Question: How much colder is it at the highest point in the area compared to the temperature at the meteostation (i.e., at 2058.1 m)?

1. About 2 degrees.

2. About 12 degrees.

3. About 7 degrees.

4. There is no difference.

Correct answers: c.

Feedback: Use aguila to animate the time series of maps with the temperature correction. Select any time step. Open the data value window by clicking on the cross and click on the map to retrieve cell values.

---

As a next step, calculate and write to disk a time series of maps containing t_i. Save the script (temp.py), run and display the new output.

We assume that precipitation falls as snow when t_iis below zero. Add statements to calculate and store a time series of maps (use the variable name freezing) that is True at a cell if the temperature is below zero degrees and False elsewhere. Use a comparison operator. Run the script and check the results.

Create two time series of maps containing respectively snowfall (m/day, variable name snowFall) and rainfall (m/day, variable name rainFall). You will need the operations ifthenelse and a Boolean operator. Write them to disk and display the results.

---

Question: Which comparison operator did you use, and where?

1. ==, in the dynamic

2. or, in the dynamic

3. while`, in the ``initial

4. > or < or >= or <=, in the dynamic

Correct answers: d.

Feedback:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('dem.map')

  def initial(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    elevationAboveMeteoStation = dem - elevationMeteoStation
    temperatureLapseRate = 0.005
    self.temperatureCorrection = elevationAboveMeteoStation * \
                                 temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

  def dynamic(self):
    precipitation = timeinputscalar('precip.tss',1)
    self.report(precipitation,'pFromTss')
    temperatureObserved = timeinputscalar('temp.tss',1)
    self.report(temperatureObserved,'tempObs')
    
    temperature = temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    self.report(freezing,'fr')
    snowFall=ifthenelse(freezing,precipitation,0.0)
    self.report(snowFall,'snF')
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)
    self.report(rainFall,'rF')

nrOfTimeSteps=181
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

---

2.4.2. The snow store

To simulate snow melt in each cell, a snow store is used containing a snow pack a_i(t) (m water equivalent) that changes over time:

a_i(t) = a_i(t-1) + s_i(t) - b_i(t)

In the equation above, (t) refers to the current time step and (t-1) to the previous time step (day). The variable s_i(t) is snowfall and b_i(t) is snowmelt. Actual snowmelt b_i(t) is modelled as:

b_i(t) = min(a_i(t), b_p,i(t))

with, b_p,i(t), potential snowmelt. Potential snowmelt is:

b_p,i(t) = mt_i, for t_i > 0

b_p,i(t) = 0, for t_i <= 0

with m, a degree day factor (here: 0.01), and t_i , the corrected termperature (see above).

Add this snow store in a stepwise manner to the model, in the following way. Let’s first assume there is no snow melt. Open temp.py and save it as snow.py. Add the variable snow representing the snow store a_i, assuming there is no snow at the start of the model run (as it is end of summer), and increasing the snowstore with snowfall in each time step.

---

Question: Without snow melt (as was assumed here), what is the maximum snowpack thickness (m water equivalent) found in the area at the end of the model run? Idem, the minimum thickness?

1. Maximum thickness: 0.11 m, minimum: 0.01 m

2. Maximum thickness: 1.12 m, minimum: 0.15 m

3. Maximum thickness: 0.85 m, minimum: 0.00 m

4. Maximum thickness: 1.96 m, minimum: 0.62 m

5. Maximum thickness: 0.57 m, minimum: 0.01 m

Correct answers: e.

Feedback:

You can get the answer by clicking on the map or by typing:

mapattr -p st000000.181

The script should look like this:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('dem.map')

  def initial(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    elevationAboveMeteoStation = dem - elevationMeteoStation
    temperatureLapseRate = 0.005
    self.temperatureCorrection = elevationAboveMeteoStation * \
                                 temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

    self.snow = 0.0

  def dynamic(self):
    precipitation = timeinputscalar('precip.tss',1)
    self.report(precipitation,'pFromTss')
    temperatureObserved = timeinputscalar('temp.tss',1)
    self.report(temperatureObserved,'tempObs')
    
    temperature= temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    self.report(freezing,'fr')
    snowFall=ifthenelse(freezing,precipitation,0.0)
    self.report(snowFall,'snF')
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)
    self.report(rainFall,'rF')

    self.snow = self.snow + snowFall
    self.report(self.snow,'snow')

nrOfTimeSteps=181
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

---

Add snowmelt to the model. First, calculate a variable potentialMelt, containing b_p,i(t). Write to disk and check the results. Second, calculate actual snowmelt b_i(t). Check the results.

---

Question: How did you calculate actual snowmelt?

1. It is equal to the potential snow melt in this case.

2. It is a fraction of the potential melt.

3. It depends on the potential melt and the actual snow thickness.

4. It is constant, so it needs to be calculated in the initial.

Correct answers: c.

Feedback:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('dem.map')

  def initial(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    elevationAboveMeteoStation = dem - elevationMeteoStation
    temperatureLapseRate = 0.005
    self.temperatureCorrection = elevationAboveMeteoStation * \
                                 temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

    self.snow = 0.0

  def dynamic(self):
    precipitation = timeinputscalar('precip.tss',1)
    self.report(precipitation,'pFromTss')
    temperatureObserved = timeinputscalar('temp.tss',1)
    self.report(temperatureObserved,'tempObs')
    
    temperature= temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    self.report(freezing,'fr')
    snowFall=ifthenelse(freezing,precipitation,0.0)
    self.report(snowFall,'snF')
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)
    self.report(rainFall,'rF')

    self.snow = self.snow + snowFall

    potentialMelt = ifthenelse(pcrnot(freezing),temperature*0.01,0)
    self.report(potentialMelt,'pmelt')
    actualMelt = min(self.snow, potentialMelt)
    self.report(actualMelt,'amelt')
    self.report(self.snow,'snow')

nrOfTimeSteps=181
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

---

And as the third step, update the snow store by subtracting actual snowmelt from the store. Write the snow store to disk (directly after subtracting actual snowmelt). Save the model under the same name snow.py.

---

Question: From what time step on is the main valley in the top-right corner snow-free, in spring? Idem, for the valley in the top-left corner?

1. Top-right corner: day 128. Top-left corner: day 130.

2. Top-right corner: day 110. Top-left corner: day 110.

3. Top-right corner: day 112. Top-left corner: day 140.

4. Top-right corner: day 92. Top-left corner: day 118.

Correct answers: a.

Feedback:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('dem.map')

  def initial(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    elevationAboveMeteoStation = dem - elevationMeteoStation
    temperatureLapseRate = 0.005
    self.temperatureCorrection = elevationAboveMeteoStation * \
                                 temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

    self.snow = 0.0

  def dynamic(self):
    precipitation = timeinputscalar('precip.tss',1)
    self.report(precipitation,'pFromTss')
    temperatureObserved = timeinputscalar('temp.tss',1)
    self.report(temperatureObserved,'tempObs')
    
    temperature= temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    self.report(freezing,'fr')
    snowFall=ifthenelse(freezing,precipitation,0.0)
    self.report(snowFall,'snF')
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)
    self.report(rainFall,'rF')

    self.snow = self.snow + snowFall

    potentialMelt = ifthenelse(pcrnot(freezing),temperature*0.01,0)
    self.report(potentialMelt,'pmelt')
    actualMelt = min(self.snow, potentialMelt)
    self.report(actualMelt,'amelt')

    self.snow = self.snow - actualMelt
    self.report(self.snow,'snow')

nrOfTimeSteps=181
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

---

2.4.3. Runoff generation

Runoff generated in each cell (m/day) is the sum of the rainfall in that cell and the snowmelt. Add a statement to snow.py calculating a time series of maps containing the total amount of runoff generated. Write results to disk and evaluate the results.

---

Question: In what time of the year is the largest amount of runoff generated?

1. Melting snow generates a large amount of runoff around time step 70.

2. Melting snow generates a large amount of runoff around time step 90.

3. Melting snow generates a large amount of runoff around time step 110.

4. Melting snow generates a large amount of runoff around time step 130.

Correct answers: d.

Feedback:

from pcraster import *
from pcraster.framework import *

class MyFirstModel(DynamicModel):
  def __init__(self):
    DynamicModel.__init__(self)
    setclone('dem.map')

  def initial(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    elevationAboveMeteoStation = dem - elevationMeteoStation
    temperatureLapseRate = 0.005
    self.temperatureCorrection = elevationAboveMeteoStation * \
                                 temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

    self.snow = 0.0

  def dynamic(self):
    precipitation = timeinputscalar('precip.tss',1)
    self.report(precipitation,'pFromTss')
    temperatureObserved = timeinputscalar('temp.tss',1)
    self.report(temperatureObserved,'tempObs')
    
    temperature= temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    self.report(freezing,'fr')
    snowFall=ifthenelse(freezing,precipitation,0.0)
    self.report(snowFall,'snF')
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)
    self.report(rainFall,'rF')

    self.snow = self.snow + snowFall

    potentialMelt = ifthenelse(pcrnot(freezing),temperature*0.01,0)
    self.report(potentialMelt,'pmelt')
    actualMelt = min(self.snow, potentialMelt)
    self.report(actualMelt,'amelt')

    self.snow = self.snow - actualMelt

    runoffGenerated = actualMelt + rainFall
    self.report(runoffGenerated,'rg')
    
    self.report(self.snow,'snow')

nrOfTimeSteps=181
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()

---