Implementing a bucket model using WorldClim layers

Species distribution models typically use climatic variables as input. The Worldclim data set includes a set of derived “Bioclim” variables that were designed to improve the interpretability of model input with respect to processes that affect organisms.
A key determinant for many plant species is length of the growing season. While this is partly determined by temperature and incoming solar radiation in temperate areas, in the tropics the growing season is limited mainly by water availability.
A conventional rule is that precipitation exceeds evapotranspiration during months with over 100mm of rainfall. This has been used as an addition to the bioclim variables when modelling. Howver this does not take into account the true water balance and introduces bias at higher latitudes and elevations.

Evapotranspiration is a complex process that integrates all climatic variables. Empirical measurements of pan evaporation and equations such as Penman Montieth produce an estimate of potential evapotranspiration (PET). Actual evapotranspiration (AET) is often much lower than this as it is affected by soil water availability. Ecophysiologists and micro-meteorologists can measure or model AET quite accurately for relatively restricted areas such as fields and some forest canopies using methods such as eddy flux covariance. However scaling up to larger spatial extents is challenging.

For the purposes of species distribution modelling precise quantative measures of AET are not necessary. The important element is to include a layer that shows relative differences in the process in a realistic manner.

One possibility would be to use Modis derived NDVI. The NDVI is correlated with productivity, as is AET, and areas with high NDVI have higher transpiraation rates than areas with low NDVI. The problem with this is that NDVI is also affected by land use. So it integrates two processes at once. There is an additional issue in irrigated areas. So while NDVI can, and should be considered in models at some point it is not the ideal input to a climate only based model.

Empirical measures of PET are available from some climate stations, but measures of pan evapotranspiration are notoriously inconsistent. This explains the absence of an evapotranspiration layer in WorldClim.

The Priestley-Taylor equation is often used to estimate PET over large areas. Although the simple form takes temperature and latitude as input, the equation does work with an estimate of net radiation so takes into account the seasonal variation of incoming radiation. It is implemented in the EcoHydRology package in R and can be adapted to run over Wordlclim layers.

## Loading required package: sp
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08 Path to GDAL shared
## files: /usr/share/gdal/1.9 Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March
## 2012, [PJ_VERSION: 470] Path to PROJ.4 shared files: (autodetected)
## OGR data source with driver: PostgreSQL 
## Source: "PG:dbname=modelling", layer: "continent"
## with 1957 features and 4 fields
## Feature type: wkbPolygon with 2 dimensions

Input data

The input to the model is a 36 band Geotiff containing the Worldclim variables minimum temperatures (months 1:12), maximum temperatures (1:12) and precipitation.

## Loading required package: raster
## raster 2.0-12 (1-September-2012)
## Loading required package: operators
## Attaching package: 'operators'
## The following object(s) are masked from 'package:base':
## options
## Loading required package: topmodel
## Loading required package: DEoptim
## DEoptim package Differential Evolution algorithm in R Authors: D. Ardia,
## K. Mullen, B. Peterson and J. Ulrich
## Loading required package: XML
clim <- brick("mintmantprec.tif")
names(clim) <- c(paste("mint", 1:12, sep = ""), paste("maxt", 1:12, sep = ""), 
    paste("prec", 1:12, sep = ""))

Calculating latitude in radians

The latitude of each pixel can be obtained from the coordinates. The function takes latitude in radians as input.

lat_rad <- coordinates(clim)[, 2] * pi/180

Calculating monthly PET

The following code takes the function in EcoHydRology and applies it to a day at the midpoint of each month (assuming 30 day months for simplicity as there is no need for false precision) .

for (i in 1:12) {
    evap <- raster(clim, 1)
    Tmax <- values(subset(clim, i + 12))/10
    Tmin <- values(subset(clim, i))/10
    d <- data.frame(day = (30 * i) - 15, Tmin, Tmax, lat_rad)
    d[] <- 0
    Es_PET <- PET_fromTemp(Jday = d$day, Tmax_C = d$Tmax, Tmin_C = d$Tmin, lat_radians = d$lat_rad) * 
    values(evap) <- Es_PET
    if (i == 1) {
        PET <<- brick(evap)
    if (i > 1) {
        PET <<- addLayer(PET, evap)

Pattern of PET

months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", 
    "Nov", "Dec")
names(PET) <- months

plot of chunk unnamed-chunk-5

Estmating AET using a simple bucket model

Actual evapotranspiration is always lower than potential evapotranspiration and can be much lower when the soil profile is well below field capacity. As soil moisture falls to close to the permanent wilting point AET will tend to zero, as no water can be extracted from the soil. The biggest challange involved in calculating AET for large areas is that soil moisture is unknown. It can never be known acurately as so many factors are involved. However soil moisture should clearly approach field capacity, and thus AET approach PET after a run of days in which rainfall exceeds evapotranspiration and fall when rain ceases. This is the basis of watershed models such as SWAT. A simple but effective way to represent this uses a “bucket model”. Once the bucket is full further rain leads to run off, and once the bucket is empty no more water can be lost. Watershed models are concerned with run off, so it is important to simulate pulses of input due to heavy rain. However a regional model can ignore this for the ppurpose at hand and simply add mean daily rainfall and subtract evapotranspiration to the bucket. AET can then be adjusted by an Alpha factor that mutiplies PET by 1 when the bucket is full and zero when empty.
In order to initialise the model a “burnin” period is needed in order to allow the bucket to find its level before monitoring the contents.

Bucket <- raster(PET, 1)

png("bucket%00d.png", height = 800, width = 900)
for (n in 1:2) {
    for (i in 1:359) {
        mn <- 1 + i%/%30
        NewAET <- raster(PET, 1)
        NewBucket <- values(Bucket)
        rain <- values(subset(clim, 24 + mn))/30
        alpha <- (NewBucket - 200)/300
        evap <- values(subset(PET, mn)) * alpha * 0.8  ##A fudge factor for stomatal control.
        NewBucket <- NewBucket + (rain) - evap
        NewBucket[NewBucket > 500] <- 500
        NewBucket[NewBucket < 200] <- 200
        values(Bucket) <- NewBucket
        values(NewAET) <- evap * (NewBucket > 200)
        if (n > 1 && (i%%30) - 15 == 0) {
            print(plot(Bucket, main = months[mn]))
            plot(cont, add = T)
            if (mn == 1) {
                AET <<- brick(NewAET)
            if (mn > 1) {
                AET <<- addLayer(AET, NewAET)
names(AET) <- months

plot of chunk unnamed-chunk-7