Bucket Hydrology Model

Monthly Hydrology for Fort Collins 2012-2013

Monthly Climate for Fort Collins 2012-2013


This is a very simple model of local surface hydrology, intended to be used for teaching and learning about weather and climate. Given monthly temperatures and precipitation for Fort Collins, CO, the model calculates potential and actual water losses by evapotranspiration and updates the level of water in a simple “bucket” representing soil moisture storage. Using the level in the bucket, the model also calculates infiltration and runoff.

The user specifies the capacity of the bucket and the initial level of stored water, both in millimeters.

Potential evapotranspiration (PET) is calculated eacxh month using the method of Thornthwaite (1948) as explained in section 11.4 of Bonan’s textbook (his equations 11.1 - 11.3):

\[ \begin{aligned} E_P = 16 \left(\frac{L}{12}\right) \left(\frac{N}{30}\right) \left(\frac{10T_a}{I}\right)^a \end{aligned} \]

where \(L\) is the average daylength during the month (hours), \(N\) is the number of days in the month, \(T_a\) is the monthly mean air temperature (Celsius),

\[ \begin{aligned} a = 6.75 x 10^{-7} I^3 - 7.71x10^{-5} I^2 + 1.79 x 10^{-2} I + 0.49 \end{aligned} \]

and \(I = \sum(\frac{T_a}{5})^{1.514}\) where the sum is carried out only over months for which \(T_a > 0\) Celsius.

Actual evapotranspiration (AET) is then calculated from PET according to

\[ \begin{aligned} AET = \beta PET \end{aligned} \]


\[ \begin{aligned} \beta = \frac{bucket(t)}{bucket_{max}} \end{aligned} \]

is the ratio of the curent level in the bucket to the bucket capacity (the amount the bucket holds when it is completely full).

In each month a fraction \(i\) of the total precipitation infiltrates into the soil bucket. We assume that this infiltration is 75% of precipitation when the bucket is empty, 25% of precipitation when the bucket is full, and that this fraction scales linearly with \(\beta\) as

The amount of water in the bucket is updated each month from

\[ \begin{aligned} bucket(t+1) = bucket(t) + (i) precip(t) - AET(t). \end{aligned} \]

Runoff is then simply given by \[ \begin{aligned} runoff(t) = (1-i) precip(t). \end{aligned} \]

In the even that the bucket “overflows,” the overcapacity amoutn is just added to runoff for that month.

The program that does the calculation is very simple, and the code that controls this website is surprisingly simple too! It is all written in the programming language R. You can read all about it on the “Website Code” tab to the right.

How this website works (including all the code!)

This website is controlled using the R package “shiny.” There are four important components:

  • The main program that updates the depth of water in the “bucket” (bucket.model.R)
  • A module that calcualtes potential evapotranspiration (thornthwaite.R)
  • A calculation of the length of the day for each month (daylength.R)
  • A program that plots the output (plot.bucket.R)
  • The user interface (ui.R) that controls the sliders and passes inputs to the model
  • The “server” (server.R) that responds to user interface events in the web browser

Scroll down or click links in the list above to read all about it!

This is the main program that updates the level in the bucket each month


bucket.model <- function(full.bucket=150, start.bucket=100) {

  # Read monthly precip and average temperature from file
  fcl <- read.table('data/FCL.monthly.txt', header=TRUE)

  # Create arrays to hold variables for later plotting
  PET <- rep(NA, 24)
  AET <- rep(NA, 24)
  runoff <- rep(NA,24)

  # Calculate potential evapotranspiration for 2012 
  # from Thornthwaite model
  PET[1:12] <- thornthwaite(fcl$temp.C[1:12], lat=41)

  # Calculate monthly potential evapotranspiration for 2013 
  # using Thornthwaite model
  PET[13:24]  <- thornthwaite(fcl$temp.C[13:24], lat=41)

  # Initialize soil moisture "bucket"
  bucket <- start.bucket  # bucket at month 0

  # Step through all the months, apportioning rainfall 
  # into runoff and infiltration, then updating the bucket
  # to calculate AET/PET in each month

  for (m in 1:24) {

    # Calculate actual evapotranspiration (AET)
    beta <- bucket / full.bucket     # soil moisture "stress"
    AET[m] <- beta * PET[m]          # Actual ET

    # Calculate infiltration fraction of precip
    #  (75% if bucket is empty, 25% when bucket is full)
    infilt <- 0.75 - 0.5 * beta

    # Runoff is what's left after infiltration into bucket
    runoff[m] <- fcl$precip.mm[m] * (1-infilt)  

    # Update the soil moisture bucket
    bucket <- bucket + fcl$precip.mm[m] - AET[m]  

    # If bucket overflows, add excess to runoff
    if (bucket > full.bucket) {
      runoff[m] <- runoff[m] + bucket - full.bucket
      bucket <- full.bucket

  # Send model output out to plotting program
  return(list(bucket=full.bucket, start=start.bucket, 
              data=data.frame(date=fcl$date, precip=fcl$precip, PET=PET, 
                              AET=AET, runoff=runoff)))


This is the module that calculates potential evapotranspiration


thornthwaite <- function(monthly.temp, lat=40.5) {

  # Use Thorthwaite (1948) relationship to calculate monthly potential evapotranspiration
  # This is Bonan's equation 11.1 and 11.2, page 161

  # input array is monthly mean temperature (Celsius), 
  #    an array of 12 values from January to December

  # second input is latitude, with a defaulyt value of 41 N

  # Only consider temperatures above freezing
  monthly.temp[monthly.temp < 0] <- 0

  # First sum a power of (positive) monthly temps to get Thornthwait's "I" value
  I <- sum((monthly.temp/5) ^ 1.514)

  # Now calculate Thornthwait's exponent "a" (Bonan Eq 11.2)
  a.power <- 6.75e-7 * I^3 - 7.71e-5 * I^2 + 1.79e-2 * I + 0.49

  # Calaulcate the average daylength (in hours) for each month
  day.Length <- rep(NA,12)
  for (month in 1:12) day.Length[month] <- daylength(lat=lat, month=month)

  # Thornthwait's "N" is the number of days in each month
  N.days <- c(31,28,31,30,31,30,31,31,30,31,30,31)

  # Finally, calculate the potential evapotranspiration (PET) for each month
  # This is Bonan's equation 11.1
  PET <- 16 * (day.Length/12) * (N.days/30) * (10 * monthly.temp / I)^a.power

  # Send the PET array back as the output of the function


Calauclate the length of the mid-month day given month and latitude


daylength <- function(lat=40., month=6){

  # Given a location and date, calculate the length of the day in hours

  # For details see http://en.wikipedia.org/wiki/Insolation


  # Orbital parameters
  ecc <- 0.0167      # Eccentricity of Earth's orbit
  eps <- 0.4091      # Obliquity (tilt) of Eart's axis (radians)
  perih <- 1.7963    # Angle between March equinox & perihelion (radians)
  equinox <- 81   # Day of the year when March equinox occurs

  # find the day of the year at mid-month
  # nDays is the number of days in each month
  nDays <- c(31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31) 
  daynum <- 0  # count days starting with zero
  if (month > 1) for (i in 1:(month-1)) daynum <- daynum + nDays[i]
  daynum <- daynum + 15

  # Theta is the position angle of the Earth in its orbit, 
  # counted in radians starting at the Vernal Equinox (March 22)
  #     theta=0 is the March Equinox, theta=90 is June Solstice
  #     theta=180 is September Equinox, theta=270 is December Solstice)
  theta <- 2 * pi * (daynum - equinox) / 365

  # Convert latitude from degrees to radians
  lat <- lat * pi / 180.

  # Declination for this time of year (theta)
  # declination is the latitude where the Sun is directly overhead @ noon
  dec <- eps * sin(theta)

  # Hour angle at sunrise & sunset
  arg <- tan(lat)*tan(dec)
  if (arg > 1) {
    h0 <- pi # perpetual summer daylight
  } else if (arg < -1) {
    h0 = 0} # perpetual winter darkness
  else {
    h0 <- acos(-arg)  # normal sunrise/sunset

  # Daylength is elapsed time between sunrise and sunset, in hours
  angular.daylength <- 2 * h0  
  daylength.hours <- 24 * angular.daylength / (2 * pi) 


This is the program that plots the output of the bucket model


plot.bucket <- function(model.output){

  # Copy model output to local arrays that extend to Jan 1 2014  
  precip <- c(model.output$data$precip, model.output$data$precip[24])
  AET <- c(model.output$data$AET, model.output$data$AET[24])
  runoff <- c(model.output$data$runoff, model.output$data$runoff[24])

  # Draw the monthly plots of precip, AET, and runoff
  plot(precip, typ='s', lwd=5, col='blue', xaxt='n', 
       ylab='mm', xlab='Month')
  lines(AET, col='red', typ='s', lwd=5)
  lines(runoff, col='black', typ='s', lwd=5)

  # Add axis labels at each mid-month
  axis(1, at=(1:25)+0.5, tick=FALSE,

  # Add tick marks to separate each month
  axis(1, at=1:25, labels=FALSE, tick=TRUE)

  # Draw a vertical line to separate the years

  # Annotate the plot
  year.pos <- 0.95 * max(precip)
  bucket.pos <- 0.9 * year.pos
  legend.x <- 10
  legend.y <- bucket.pos

  text(7, year.pos, '2012', cex=1.5)
  text(19, year.pos, '2013', cex=1.5)
  text(4, bucket.pos, cex=1.1, paste('Bucket size =',model.output$bucket))
  text(4, 0.9*bucket.pos, cex=1.1, paste('Initial bucket =',model.output$start))

  legend(x=legend.x, y=legend.y, c('Precip', 'AET', 'runoff'),
         col=c('blue','red','black'), lwd=c(5,5,5))

This is the user interface that actually controls this website



# Define UI for the bucket model of surface hydrology

  #  Application title
  headerPanel("Bucket Hydrology Model"),

  # Sidebar with sliders that control various available options

    sliderInput("full.bucket", "Bucket Depth:", 
                min = 10, max = 200, value = 150, step=10),

    sliderInput("start.bucket", "Initial Level:", 
                min = 0, max = 200, value = 125, step=5)


   # Main panel consists of a tabset displays model output or model description
                h3("Monthly Hydrology for Fort Collins 2012-2013"),
                tags$style(type="text/css", ".tab-content { overflow: visible; }")), 
       tabPanel("Input Data",
                h3("Monthly Climate for Fort Collins 2012-2013"),
       tabPanel("Website Code",

Here is the server module that responds to user interface events



# Source required R scripts

shinyServer(function(input, output) {

  # Monthly hydrology plot
  output$hydrology.plot <- renderPlot(

  # Table of input climate data
  output$input.data <- renderTable(
    read.table('data/FCL.monthly.txt', header=TRUE), include.rownames=FALSE)