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} \]
where
\[ \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.
This website is controlled using the R package “shiny.” There are four important components:
Scroll down or click links in the list above to read all about it!
bucket.model.R:
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)))
}
thornthwaite.R:
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
return(PET)
}
daylength.R:
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)
return(daylength.hours)
}
plot.bucket.R:
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,
labels=c(rep(c('J','F','M','A','M','J','J','A','S','O','N','D'),2),''))
# Add tick marks to separate each month
axis(1, at=1:25, labels=FALSE, tick=TRUE)
# Draw a vertical line to separate the years
abline(v=13)
# 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))
}
ui.R:
library(shiny)
library(markdown)
# Define UI for the bucket model of surface hydrology
shinyUI(pageWithSidebar(
# Application title
headerPanel("Bucket Hydrology Model"),
# Sidebar with sliders that control various available options
sidebarPanel(
img(src='bucket.jpg'),
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
mainPanel(
tabsetPanel(
tabPanel("Output",
h3("Monthly Hydrology for Fort Collins 2012-2013"),
plotOutput("hydrology.plot"),
tags$style(type="text/css", ".tab-content { overflow: visible; }")),
tabPanel("Input Data",
h3("Monthly Climate for Fort Collins 2012-2013"),
tableOutput("input.data")),
tabPanel("Description",
includeHTML('doc/bucket.description.html')),
tabPanel("Website Code",
includeMarkdown('doc/bucket.code.md'))
)
)
))
server.R:
library(shiny)
# Source required R scripts
source('model/daylength.R')
source('model/thornthwaite.R')
source('model/bucket.model.R')
source('model/plot.bucket.R')
shinyServer(function(input, output) {
# Monthly hydrology plot
output$hydrology.plot <- renderPlot(
plot.bucket(bucket.model(full.bucket=input$full.bucket,
start.bucket=input$start.bucket))
)
# Table of input climate data
output$input.data <- renderTable(
read.table('data/FCL.monthly.txt', header=TRUE), include.rownames=FALSE)
})