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.


How this website works (including all the code!)


  • 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)


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)

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


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

  # Declination for this time of year (theta)
  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)


                h3("Monthly Hydrology for Fort Collins 2012-2013"),
       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)