Modelling Framework

Introduction

The model development cycle consists typically of two steps: the conversion of the conceptual model structure into a numerical implementation, and the assimilation or calibration of the model with observed data (Karssenberg, in prep). Most of the available tools are either dedicated to model development (e.g. ModelMaker, Stella, Matlab), or to calibration or assimilation (e.g. PEST, ReBEL). The desired combination of those tools usually requires an additional software layer.

This framework offers a combined interface for the task of model construction and optimisation. Spatial operations can be used from the PCRaster library - temporal, deterministic and stochastic modelling capabilities are added with this framework. Modelling options are Monte Carlo simulations (e.g. Doucet et al, 2001), Particle filter (e.g. Xiong, 2006) and Ensemble Kalman filter (e.g Simon, 2006).

Example

The example script shows a spatio-temporal model executed in a Monte Carlo simulation. The process descriptions simulating the snow pack of a catchment are given in the initial and dynamic sections of the SnowModel class. The Monte Carlo framework executes the premcloop to generate initial values valid for all samples, and afterwards 100 realisations each running 180 time steps. Finally the postmcloop calculates statistics of the ensemble results.

For a more detailed description of the model see Karssenberg et al (2010) and the manual.

from pcraster import *

class SnowModel(object):
  def __init__(self):
    setclone("clone.map")

  def premcloop(self):
    dem = self.readmap("dem")
    self.ldd = lddcreate(dem, 1e31, 1e31, 1e31, 1e31)
    elevationMeteoStation = scalar(2058.1)
    self.elevationAboveMeteoStation = dem - elevationMeteoStation
    self.degreeDayFactor = 0.01

  def initial(self):
    self.snow = scalar(0)
    self.temperatureLapseRate = 0.005 + (mapnormal() * 0.001)
    self.report(self.temperatureLapseRate, "lapse")
    self.temperatureCorrection = self.elevationAboveMeteoStation * self.temperatureLapseRate

  def dynamic(self):
    temperatureObserved = self.readDeterministic("tavgo")
    precipitationObserved = self.readDeterministic("pr")
    precipitation = max(0, precipitationObserved * (mapnormal() * 0.2 + 1.0))
    temperature = temperatureObserved - self.temperatureCorrection
    snowFall = ifthenelse(temperature < 0, precipitation, 0)
    self.snow = self.snow + snowFall
    potentialMelt = ifthenelse(temperature > 0, temperature * self.degreeDayFactor, 0)
    actualMelt = min(self.snow, potentialMelt)
    self.snow = max(0, self.snow - actualMelt)
    rain = ifthenelse(temperature >= 0, precipitation, 0)
    discharge = accuflux(self.ldd, actualMelt + rain)
    self.report(self.snow, "s")
    self.report(discharge, "q")

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

myModel = snowModel()
dynamicModel = DynamicFramework(myModel, 180)
mcModel = MonteCarloFramework(dynamicModel, 100)
mcModel.run()

Simulation results

The image below shows a comparison of model results simulated in Monte Carlo mode (upper row) and with the particle filter (lower row, filter moments at time steps 70, 100 and 150). The left images show the snow pack (m) results of 10 realisations, the centre image the map view with currently selected cell and the right images show the cumulative probability density function.

Distribution of software

The software is part of the current release of PCRaster.

Acknowledgements

This research was supported by the Dutch research programme Space for Geo-Information (Ruimte voor Geo-Informatie, RGI) as RGI project 313 ("Online coupling of spatial optimisation tools and spatially distributed simulation models") with the cooperation partners PCRaster Environmental Software, Department of Physical Geography, Utrecht University and the Joint Research Centre.

References

Doucet, A., de Freitas, N. and Gordon, N. (2001): Sequential Monte Carlo Methods in Practice. Springer Xiong, X., Navon, I. M. and Uzunoglu, B. (2006): A note on the particle filter with posterior Gaussian resampling. Tellus Series A: Dynamic Meteorology and Oceanography Simon, D. (2006): Optimal State Estimation: Kalman, H and Nonlinear Approaches. Wiley Interscience

Karssenberg, D., Schmitz, O., Salamon, P., De Jong, K. and Bierkens, M.F.P., 2010, A software framework for construction of process-based stochastic spatio-temporal models and data assimilation. Environmental Modelling & Software, 25(4), pp. 489-502, Link.