# Bucket Hydrology Model

### MODEL DESCRIPTION

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.

### 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

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

bucket.model.R:

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

# Read monthly precip and average temperature from file

# 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.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)

}


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

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


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

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


#### This is the user interface that actually controls this website

ui.R:

library(shiny)
library(markdown)

# Define UI for the bucket model of surface hydrology
shinyUI(pageWithSidebar(

#  Application title

# 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')) ) ) ))  #### Here is the server module that responds to user interface events 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(