Water Balance Calculations

WaterBalance.R

# calculating water balance based on Willmott, Rowe, and Mintz 1985 J. Clim.

Temp <- read.csv(file.choose()) # monthly_temp.csv
Precip <- read.csv(file.choose()) # monthly_precip.csv
Temp <- Temp[Temp[,1]>1905,]
Precip <- Precip[Precip[,1]>1905,]

T <- t(Temp)[-1,]
colnames(T) <- Temp[,1]
for(r in 1:12) T[r,is.na(T[r,])] <- round(mean(T[r,!is.na(T[r,])]),2)
P <- t(Precip)[-1,]
colnames(P) <- Precip[,1]
P[is.na(P)] <- 0
# check matching years
if(any(Temp[,1] != Precip[,1])) print("ERROR: non-matching years")

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(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)
nd[2,which(Temp[,1]%in%leapyears)] <- 29

nh <- c(9.75,10.75,11.95,13,25,14.32,14.87,14.6,13.67,12.42,11.18,10.05,9.47) # number of daylight hours on the 15th of the month
wstar <- 300 # soil water storage capacity, mm
w0 <- ws0 <- 0 # initial water... set 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)
 } 
}


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=T)

plot(PET[2,],type="o",col="royalblue",ylim=c(0,max(PET)),xlab="Month",ylab="Moisture (mm)")
points(AET[2,]~seq(0.95,11.95,1),type="o",col="deeppink")
points(PPT[2,]~seq(1.05,12.05,1),type="o",col="darkgreen")
for(i in 1:12){
 lines(c(i,i),c(PET[1,i],PET[3,i]),col="royalblue",lty=3)
 lines(c(i-0.05,i-0.05),c(AET[1,i],AET[3,i]),col="deeppink",lty=3)
 lines(c(i+0.05,i+0.05),c(PPT[1,i],PPT[3,i]),col="darkgreen",lty=3)
}

legend("topright",c("PET","AET","PPT"),lty=1,col=c("royalblue","deeppink","darkgreen"),bty="n")

WB <- cbind(Temp[,1],apply(Eo-Ea,2,sum),apply(Ea,2,sum))
colnames(WB) <- c("Year","Deficit","AET")