3.7. Dynamic stochastic modelling: snow melt model

3.7.1. Snow melt as a stochastic model

In the dynamic modelling exercises you created a deterministic snow melt model. Here we will built upon this work by propagating input errors to its outputs. As a start, the data set contains a script that is exactly the same as the script created in the deterministic modelling exercises, apart from that the script has been modified to be able to do stochastic modelling. Also, we will use smaller input maps to speed up calculations.

Go to the folder snow in the folder containing your exercise data, and open the script snow.py. Be sure to remind what calculations were included (if needed, have a look at the deterministic modelling exercises). It uses the stochastic dynamic modelling framework, and thus it will run for 50 Monte Carlo loops (as you can see at the bottom). However in the snow.py, no stochastic inputs are used (yet).

Run the script. Display the input timeseries, two realizations of snow cover and discharge, and check whether they are the same. Also display dem.map and you will see that we are using a smaller area now.

3.7.2. Stochastic dynamic input: precipitation

To represent uncertainty in precipitation, we define the following error model:

p(t) = pm(t) x e(t)

With pm(t) the observed precipitation, p(t) the stochastic variable precipitation (including uncertainty), and e(t) a random non-spatial (i.e., the same for all cells) variable for each time step t with mean one and variance 0.04.

Open snow.py and save as precip.py. Modify precip.py to represent rainfall as a stochastic variable. You need to use the function mapnormal. And note that you will need to convert variance to standard deviation.

Save the script and run the model. Display 12 realizations of precipitation and snow cover. You need the scenarios, multi, and timesteps Aguila options. Animate them, and create timeseries for different locations.

Question: Do you think our error model results in a large error in the total amount of rain over the whole modelling period? Explain.

  1. No, the errors are assigned each time step thus resulting in a large error (additive error).

  2. Yes, the errors are assigned each time step thus resulting in a large error (additive error).

  3. No, the errors cancel each other out as they are assigned each time step independently.

  4. Yes, the errors cancel each other out as they are assigned each time step independently.

Correct answers: c.

Feedback: None

To be able to create plots of probability density distributions, you can calculate the mean, variance, and percentiles in the postmcloop. Do this for precipitation, snow depth (snow), and discharge (q). To reduce the calculation time, use a limited set of percentiles:

percentiles=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]

Save the precip.py script and run. Note that it will take some time at the end to calculate the percentiles. If the script is really too slow, you can consider using less Monte Carlo samples (e.g. 25 instead of 50), but note that this will result in a less precise calculation of the errors.

Use aguila to display and animate probability density distributions of precipitation, snow depth, and discharge, together with the map locations.map. Also, plot the variance of these attributes using Aguila, together with locations.map.

Question: What is the variance in snow depth at the two locations on locations.map, for time step 120 and 170 days? Write it down too.

  1. Day 120: Location 1, 6.16, Location 2, 1.55; Day 170: Location 1, 7.90; Location 2, 1.12

  2. Day 120: Location 1, 8.06 x 10-5, Location 2, 0.36 x 10-7; Day 170: Location 1, 0.12; Location 2, 0.0

  3. Day 120: Location 1, 4.99, Location 2, 1.41; Day 170: Location 1, 6.44; Location 2, 1.35

  4. Day 120: Location 1, 3.06 x 10-5, Location 2, 2.32 x 10-7; Day 170: Location 1, 0.0; Location 2, 0.0

Correct answers: d.

Feedback: None

Question: What is the variance in discharge at the two locations on locations.map, for time step 120 and 170 days? Write it down too.

  1. Day 120: Location 1, 0.01, Location 2, 0.00; Day 170: Location 1, 1.75 * 10^10; Location 2, 3.87 * 10^13

  2. Day 120: Location 1, 0.04, Location 2, 0.00; Day 170: Location 1, 2.75 * 10^4; Location 2, 1.61 * 10^9

  3. Day 120: Location 1, 0.09, Location 2, 0.00; Day 170: Location 1, 3.75 * 10^10; Location 2, 9.61 * 10^13

  4. Day 120: Location 1, 0.00, Location 2, 0.00; Day 170: Location 1, 7.83 * 10^4; Location 2, 2.52 * 10^9

Correct answers: d.

Feedback:

from pcraster import *
from pcraster.framework import *

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

  def premcloop(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    self.elevationAboveMeteoStation = dem - elevationMeteoStation
    self.ldd=lddcreate(dem,1e31,1e31,1e31,1e31)
    self.report(self.ldd,'ldd')

  def initial(self):
    temperatureLapseRate = 0.005
    self.temperatureCorrection = self.elevationAboveMeteoStation * temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

    self.snow=0.0

  def dynamic(self):
    #precipitation = timeinputscalar('precip.tss',1)
    precipitationFactor=max(0, mapnormal()*0.2 + 1.0)
    precipitation = timeinputscalar('precip.tss',1) * precipitationFactor
    
    self.report(precipitation,'p')
    temperatureObserved = timeinputscalar('temp.tss',1)
    temperature= temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    snowFall=ifthenelse(freezing,precipitation,0.0)
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)

    self.snow = self.snow+snowFall

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

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

    runoffGenerated = actualMelt + rainFall

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

  def postmcloop(self):
    names=['p','q', 'snow']
    sampleNumbers=self.sampleNumbers()
    timeSteps=self.timeSteps()
    percentiles=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
    mcaveragevariance(names,sampleNumbers,timeSteps)
    mcpercentiles(names,percentiles,sampleNumbers,timeSteps)

nrOfTimeSteps=181
nrOfSamples=50
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
mcModel = MonteCarloFramework(dynamicModel, nrOfSamples)
# dynamicModel.run()
mcModel.run()

3.7.3. Stochastic parameters: temperature lapse rate

To represent uncertainty in the temperature lapse rate parameter, we replace the deterministic value of l by a non-spatial (i.e. the same for all grid cells) stochastic variable with mean 0.005 and variance 0.000001.

Open precip.py and save as lapse.py. Modify the script using the stochastic l instead of the deterministic fixed value. As the uncertainty in lapse rate has effect on the uncertainty in temperature, also calculate the mean, variance, and percentiles for temp. Run the script.

Display the probability density function of temperature, to see the effect of uncertain lapse rate on temperature.

Question: Write down the variance in snow depth at the two locations on locations.map, for time step 120 and 170 days. Compare the result with the version of the model that ignored uncertainty in lapse rate.

  1. Day 120: Location 1, 3.91 x 10-5, Location 2, 8.02 x 10-7; Day 170: Location 1, 0; Location 2, 0

  2. Day 120: Location 1, 4.99, Location 2, 1.41; Day 170: Location 1, 6.44; Location 2, 0.35

  3. Day 120: Location 1, 3.94 x 10-8, Location 2, 6.02 x 10-8; Day 170: Location 1, 0; Location 2, 0

  4. Day 120: Location 1, 5.16, Location 2, 4.41; Day 170: Location 1, 1.75; Location 2, 6.28

Correct answers: a.

Feedback:

from pcraster import *
from pcraster.framework import *

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

  def premcloop(self):
    dem = self.readmap('dem')
    elevationMeteoStation = 2058.1
    self.elevationAboveMeteoStation = dem - elevationMeteoStation
    self.ldd=lddcreate(dem,1e31,1e31,1e31,1e31)
    self.report(self.ldd,'ldd')

  def initial(self):
    #temperatureLapseRate = 0.005
    temperatureLapseRate = 0.005 + mapnormal()*0.001
    self.temperatureCorrection = self.elevationAboveMeteoStation * temperatureLapseRate
    self.report(self.temperatureCorrection,'tempCor')

    self.snow=0.0

  def dynamic(self):
    precipitationFactor=max(0, mapnormal()*0.2 + 1.0)
    precipitation = timeinputscalar('precip.tss',1) * precipitationFactor
    self.report(precipitation,'p')

    temperatureObserved = timeinputscalar('temp.tss',1)
    temperature= temperatureObserved - self.temperatureCorrection
    self.report(temperature,'temp')

    freezing=temperature < 0.0
    snowFall=ifthenelse(freezing,precipitation,0.0)
    rainFall=ifthenelse(pcrnot(freezing),precipitation,0.0)

    self.snow = self.snow+snowFall

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

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

    runoffGenerated = actualMelt + rainFall

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

  def postmcloop(self):
    names=['p','q', 'snow','temp']
    sampleNumbers=self.sampleNumbers()
    timeSteps=self.timeSteps()
    percentiles=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
    mcaveragevariance(names,sampleNumbers,timeSteps)
    mcpercentiles(names,percentiles,sampleNumbers,timeSteps)

nrOfTimeSteps=181
nrOfSamples=50
myModel = MyFirstModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
mcModel = MonteCarloFramework(dynamicModel, nrOfSamples)
mcModel.run()