2.5. Neighbourhood operations with defined topology: the snow melt model

2.5.1. The local drain direction map

In the previous exercise you created a snowmelt model. Your model calculated for each time step (day) the total amount of runoff that was generated in a cell. Here, you will transport this runoff downstream. One approach is to assume that runoff occurs towards a neighbouring cell with the steepest downhill path. These directions are stored in a so-called local drain direction map.

Open snow.py and save it as runoff.py. Add a lddcreate statement calculating the local drain direction map. Run the model and have a look at the local drain direction map.

---

Question: Where in the script did you enter the lddcreate function? Explain.

1. At the top.

2. In the __init__.

3. In the initial.

4. In the dynamic

Correct answers: c.

Feedback:

You need the following lines in the initial. The lddcreate function is added to the initial section, because it needs to be calculated only once, at the start. It needs to be defined as a global variable to be later used in the dynamic section for routing of water, so this is why a variable name starting with self. is used.

    self.ldd=lddcreate(dem,1e31,1e31,1e31,1e31)
    self.report(self.ldd,'ldd')

---

2.5.2. Drain all runoff within one time step: the accuflux function

When the size of a study area is sufficiently small and the time step duration sufficiently long that it can be assumed that all runoff generated in a time step reaches the outflow point in the same time step, runoff can be transported downstream using an accu function.

Add an accuflux statement to runoff.py calculating discharge for each grid cell. You need to calculate discharge as a result of the sum of rainfall (m/day) and snowmelt (m/day). Discharge should be calculated as m3/day. You can get the area of a cell (here in m3 as the length of the cells on the map is entered in metres when creating the map) with the function cellarea. Save the runoff.py script and run.

Animate the discharge using Aguila. Try a logarithmic scale for better visualisation (right-click on the legend, select ‘Edit draw properties’, set ‘Colour assignment’ to logarithmic and set the ‘minimum cutoff’ to a positive (non zero) value).

---

Question: What is the maximum discharge in the study area?

1. 0.3·10^6 m3/day

2. 1.3·10^6 m3/day

3. 2.3·10^6 m3/day

4. 1.3·10^7 m3/day

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

    self.ldd=lddcreate(dem,1e31,1e31,1e31,1e31)
    self.report(self.ldd,'ldd')

  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,'freez')
    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')

    discharge=accuflux(self.ldd,runoffGenerated*cellarea())
    self.report(discharge,'q')

    self.report(self.snow,'snow')

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

---