Code for calculating annual water balance with data from NCDC

BerdanierKwit_waterbalance.R

# Calculating water balance based on algorithms in Willmott, Rowe, and Mintz (1985) J. Clim. 5:589-606
# Aaron Berdanier and Matt Kwit (2012)
# Questions or comments to: aaron.berdanier@duke.edu

# Packages:
library(reshape)

# Load the function:
waterbalance <- function(rawdat,nh,wstar){

 redat <- rawdat[rawdat$STATION_NAME==stations[choice],]
 yearmo <- cbind(as.numeric(substr(redat$DATE,1,4)),as.numeric(substr(redat$DATE,5,6)))
 redat <- redat[yearmo[,1]%in%min(yearmo[yearmo[,2]==1,1]):max(yearmo[yearmo[,2]==12,1]),]
 yearmo <- yearmo[yearmo[,1]%in%min(yearmo[yearmo[,2]==1,1]):max(yearmo[yearmo[,2]==12,1]),]

 meltd <- cbind(yearmo,redat[,4:5])
 colnames(meltd) <- c("Year","Mo","P","T")
 T <- as.matrix(cast(meltd,Mo~Year,function(x) mean(x)/10,value="T"))
 P <- as.matrix(cast(meltd,Mo~Year,function(x) mean(x)/10,value="P"))
 T[is.nan(T)] <- NA
 P[is.nan(P)] <- NA

 # T <- as.matrix(T)
 # P <- as.matrix(P)

 # Missing values are given the long-term average
 # could be changed for missing values to interpolate surrounding values too... might be better
 for(r in 1:12) T[r,is.na(T[r,])] <- round(mean(T[r,!is.na(T[r,])]),2)
 # Missing values assigned zero
 P[is.na(P)] <- 0

 ny <- ncol(P)	# number of observation years
 Eo <- Eox <- Ea <- w15 <- M15 <- matrix(NA,12,ny)

 nd <- matrix(rep(c(31,28,31,30,31,30,31,31,30,31,30,31),ny),12,ny,byrow=F) 	# number of days in each month, year
 leapyears <- c(1900,1904,1908,1912,1916,1920,1924,1928,1932,1936,1940,1944,1948,1952,1956,1960,1964,1968,1972,1976,1980,1984,1988,1992,1996,2000,2004,2008,2012)
 nd[2,which(as.numeric(colnames(T))%in%leapyears)] <- 29

 w0 <- ws0 <- 0 # set initial water storage to zero

 for(y in 1:ny){
  I <- sum(!is.nan((T[,y]/5)^1.514))
  a <- 6.75E-7*I^3 - 7.71E-5*I^2 + 1.79E-2*I + 0.49

  for (m in 1:12){

   # Potential evapotranspiration  
   if(T[m,y] < 0){
    Eox[m,y] <- 0
   }else{
    if(T[m,y] < 26.5){
     Eox[m,y] <- 16*(10*T[m,y]/I)^a
    }else{
     Eox[m,y] <- -415.85 + 32.24*T[m,y] - 0.43*T[m,y]^2
    }
   }
   Eo[m,y] <- Eox[m,y]*((nd[m,y]/30)*(nh[m]/12))

   D <- B <- S <- M <- numeric(nd[m,y])
   Pd <- P[m,y]/nd[m,y]
   Ps <- Pr <- 0

   w <- ws <- numeric(nd[m,y]+1)
   w[1] <- w0
   ws[1] <- ws0

   for (d in 1:nd[m,y]){
    if(T[m,y] < -1){
     Ps <- Pd
    }else{
     Pr <- Pd
    }
   
    # Snowmelt per day
    M[d] <- 2.63 + 2.55*T[m,y] + 0.0912*T[m,y]*Pr
    if(M[d] < 0) M[d] = 0 # snowmelt cannot be negative
    if(M[d] > ws[d] + Ps) M[d] = ws[d] + Ps # snowmelt cannot be greater than available water
    ws[d+1] <- ws[d] + Ps - M[d] # SWE for next day

    # Demand per day
    D[d] <- M[d] + Pd - Eo[m,y]/30
    if(D[d] < 0){ # evaporative demand
     B[d] <- 1-exp(-6.68*w[d]/wstar)
    }else{ # recharge
     B[d] <- 1
    }

    # soil water at the end of the day
    w[d+1] <- w[d] + B[d]*D[d]
    if(w[d+1] > wstar){ # surplus!
     S[d] <- w[d+1] - wstar
     w[d+1] <- wstar
    }
   }

   # Change in soil water over the month
   dw <- tail(w,n=1) - w[1]
   w15[m,y] <- w[16] # soil moisture on the 15th of the month
   M15[m,y] <- M[15]

   Ea[m,y] <- P[m,y] + sum(M) - dw - sum(S) 

   w0 <- tail(w,n=1) # assign new w0
   ws0 <- tail(ws,n=1)
  } 
 }
 Eo <<- Eo
 Ea <<- Ea
 P <<- P
}

# Steps:
# 1. Download data from http://www.ncdc.noaa.gov/cdo-web/#t=secondTabLink
#  - Choose "Monthly Summaries GHCND"
#  - Identify stations and download temperature and precipitation data
# 2. Load csv...
rawdat <- read.csv(file.choose())
head(rawdat)
rawdat[rawdat==9999] <- NA
stations <- unique(rawdat$STATION_NAME)
print("Choose a station...");matrix(stations)
# 3. Choose a station.
choice <- 3

# 4. Need number of daylight hours on the 15th of each month at the site
# http://aa.usno.navy.mil/data/docs/Dur_OneYear.php
nh <- 
# Fort Collins, CO
 c(9.6,10.66666667,11.95,13.31666667,14.46666667,15.06666667,14.78333333,13.76666667,12.45,11.11666667,9.916666667,9.283333333)
# Manhattan, KS
# c(9.716666667,10.73333333,11.95,13.25,14.35,14.91666667,14.65,13.68333333,12.43333333,11.16666667,10.01666667,9.433333333)
# Minot ND
# c(8.783333333,10.23333333,11.9,13.7,15.23333333,16.06666667,15.66666667,14.3,12.56666667,10.83333333,9.216666667,8.366666667)

# 5. Need soil water storage capacity, mm
wstar <- 300

waterbalance(rawdat,nh,wstar)

PET <- apply(Eo,1,quantile,c(0.1,0.5,0.9))
AET <- apply(Ea,1,quantile,c(0.1,0.5,0.9))
PPT <- apply(P,1,quantile,c(0.1,0.5,0.9),na.rm=TRUE)