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 Celsius) to the map variable temperatureObserved. Write the results to disk, run the model, and have a look at the output.

---

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

1. open

2. timeinputscalar

3. timeinputnominal

4. timeinput

5. report

Correct answer: 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 12: 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 answer: 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.

import pcraster as pcr
import pcraster.framework as pcrfw


class MyFirstModel(pcrfw.DynamicModel):
    def __init__(self):
        pcrfw.DynamicModel.__init__(self)
        pcr.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 = pcr.timeinputscalar("precip.tss", 1)
        self.report(precipitation, "pFromTss")
        temperatureObserved = pcr.timeinputscalar("temp.tss", 1)
        self.report(temperatureObserved, "tempObs")


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

---

Question 13: 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 answer: 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 14: 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 answer: d

Feedback: None

import pcraster as pcr
import pcraster.framework as pcrfw


class MyFirstModel(pcrfw.DynamicModel):
    def __init__(self):
        pcrfw.DynamicModel.__init__(self)
        pcr.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 = pcr.timeinputscalar("precip.tss", 1)
        self.report(precipitation, "pFromTss")
        temperatureObserved = pcr.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 = pcr.ifthenelse(freezing, precipitation, 0.0)
        self.report(snowFall, "snF")
        rainFall = pcr.ifthenelse(pcr.pcrnot(freezing), precipitation, 0.0)
        self.report(rainFall, "rF")


nrOfTimeSteps = 181
myModel = MyFirstModel()
dynamicModel = pcrfw.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 15: 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 answer: e

Feedback: You can get the answer by clicking on the map or by typing: mapattr -p st000000.181

import pcraster as pcr
import pcraster.framework as pcrfw


class MyFirstModel(pcrfw.DynamicModel):
    def __init__(self):
        pcrfw.DynamicModel.__init__(self)
        pcr.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 = pcr.timeinputscalar("precip.tss", 1)
        self.report(precipitation, "pFromTss")
        temperatureObserved = pcr.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 = pcr.ifthenelse(freezing, precipitation, 0.0)
        self.report(snowFall, "snF")
        rainFall = pcr.ifthenelse(pcr.pcrnot(freezing), precipitation, 0.0)
        self.report(rainFall, "rF")

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


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

---

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 16: 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 answer: c

Feedback: None

import pcraster as pcr
import pcraster.framework as pcrfw


class MyFirstModel(pcrfw.DynamicModel):
    def __init__(self):
        pcrfw.DynamicModel.__init__(self)
        pcr.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 = pcr.timeinputscalar("precip.tss", 1)
        self.report(precipitation, "pFromTss")
        temperatureObserved = pcr.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 = pcr.ifthenelse(freezing, precipitation, 0.0)
        self.report(snowFall, "snF")
        rainFall = pcr.ifthenelse(pcr.pcrnot(freezing), precipitation, 0.0)
        self.report(rainFall, "rF")

        self.snow = self.snow + snowFall

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


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

---

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 17: 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 answer: a

Feedback: None

import pcraster as pcr
import pcraster.framework as pcrfw


class MyFirstModel(pcrfw.DynamicModel):
    def __init__(self):
        pcrfw.DynamicModel.__init__(self)
        pcr.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 = pcr.timeinputscalar("precip.tss", 1)
        self.report(precipitation, "pFromTss")
        temperatureObserved = pcr.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 = pcr.ifthenelse(freezing, precipitation, 0.0)
        self.report(snowFall, "snF")
        rainFall = pcr.ifthenelse(pcr.pcrnot(freezing), precipitation, 0.0)
        self.report(rainFall, "rF")

        self.snow = self.snow + snowFall

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

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


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

---

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 18: 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 answer: d

Feedback: None

import pcraster as pcr
import pcraster.framework as pcrfw


class MyFirstModel(pcrfw.DynamicModel):
    def __init__(self):
        pcrfw.DynamicModel.__init__(self)
        pcr.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 = pcr.timeinputscalar("precip.tss", 1)
        self.report(precipitation, "pFromTss")
        temperatureObserved = pcr.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 = pcr.ifthenelse(freezing, precipitation, 0.0)
        self.report(snowFall, "snF")
        rainFall = pcr.ifthenelse(pcr.pcrnot(freezing), precipitation, 0.0)
        self.report(rainFall, "rF")

        self.snow = self.snow + snowFall

        potentialMelt = pcr.ifthenelse(pcr.pcrnot(freezing), temperature * 0.01, 0)
        self.report(potentialMelt, "pmelt")
        actualMelt = pcr.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 = pcrfw.DynamicFramework(myModel, nrOfTimeSteps)
dynamicModel.run()

---