This is a report on how to calculate actual evapotranspiration of Rock Creek Meadow located in the Eastern Sierra’s of California. There are three methods we will use:

The Priestly-Taylor evapotranspiration equation is listed as:

\[ET = α \frac{( ∆ \ (Rn-G))} {(λ_v\ (\ ∆+γ \ ))}\] \[ α = \ Aerodynamic \ canopy \ resistance\] \[ ∆ = \ Slope \ of \ saturated \ vapor \ pressure \ temperature \ curve \ (\frac{kPa}{c})\] \[ λ_v = \ Volumemetric \ latent\ heat\ of \ vaporization (\frac{2453\ MJ}{m^{3}})\] \[ R_n =\ Net \ Radiation\ (\frac{MJ}{m^{2}}{day^{-1}})\] \[ G = \ Soil \ Heat \ Flux\ Density \ (\frac{MJ}{m^{2}}{day^{-1}})\] \[ γ = \ Psychometric \ constant (\frac{kPa}{C})\]

Background

Priestly - Taylor Method

This equation was developed as a substitute for the commonly used Pennman Monteith equation by Priestly and Taylor in 1972. It is a simplified version of the Pennman Monteith equation, that requires less observational data input, as only solar radiation and temperature inputs are required. This allows for the addition of an empirically derived constant, α, which is vapor deficit and convection terms in a single constant.

In the article “Estimation and comparison from MODIS and AVHRR sensors for clear sky days over the Southern Great Plains”, Namrata Batra and colleagues discuss remote sensing techniques and MODIS NDVI data to determine actual evapotranspiration (Batra et al. 2006).

NDVI data was used to linearly interpolated α (aerodynamic canopy resistance) values. They found a range of values for aerodynamic canopy resistance variable, as there was a relationship between NDVI, surface temperature and aerodynamic canopy resistance that was discovered (Jiang and Islam 2001). s

To find this relationship, they used multiple remotely sensed images and linear interpolation to find variables within the Priestly-Taylor equation stating, “we begin with the detailed analysis of remotely sensed parameters: To and NDVI, followed by the coefficient ‘Φ’ which gives us an estimate of EF” (note that Φ is the same variable as α).

This methodology allowed them to determine AET over a large-scale area using remotely sensed imagery. Batra’s paper served to provide the basis of applying NDVI data to our study area in Rock Creek and a similar methodology of Batra to obtain AET values.

One limitation of our data is that is is missing soil heat flux density values, meaning we had to find a method for determining soil heat flux without the actual measured values. The solution was to use the assumption that about 10 percent of solar radiation recorded each day equates to soil heat flux levels, or utilize Batra’s method which uses solar radiation and NDVI data to calculate heat flux.

Another limitation is that NDVI data is compiled over 16-day average, but the climate and atmospheric data occurs every 30 minutes of each day. The decision was made to compile our climate and atmospheric data into 16 days intervals, so calculations with the NDVI will be consistent to the dates of the climate and atmospheric data.

Data Sources and More Information:

The climate, atmospheric, and additional data on Rock Creek data was provided by Chris Surfleet, a hydrology professor in the Natural Resources and Environmental Sciences department at Cal Poly San Luis Obispo.

To learn more about the ET equation, please visit Jason Mercer’s discussion:
Penman Monteith and Priestley Taylor ET

The NDVI data was found using APPEARS. Using a shapefile of Rock Creek, we extracted MODIS 16-day 200 meter NDVI data which was converted into a statistical csv.

Data Manipulation and Variable Calculations

Two csv files of atmospheric data including absolute pressure and temperature were provided for June - September of 2019 and then September 2019 - July 2020. For ease of compilation, these atmospheric csv files were combined into one csv file which was read into RStudio. This allowed for easy readability and less time spent managing dates and times within RStudio.

A c;imate csv data was also provided, which included date time, solar radiation metrics, wind speed (mph), gust speed, temperature relative humidity, rain and snow depth, and current. Data collection occurs every 30 minutes each day at Rock Creek Climate station, providing accurate and consistent measurements for the given time period. Before reading in the files, some pre-processing and clean up was required to ensure correct implementation into RStudio.

The packages used are:

library(ggplot2)
library(dplyr)
library(lubridate)
library(readr)
library(readxl)
library(tidyverse)
library(cowplot)
library(rmarkdown)
library(knitr)
library(Evapotranspiration)
library(zoo)
library(tidyr)

Let’s read in the climate, NDVI, and atmospheric data now:

climateRC   <-  read_excel("~/Documents/NR339/Climate/excel/climateRC.xls",  
                        col_types = c("numeric", "date", "numeric", 
                                      "numeric", "numeric", "numeric", 
                                      "numeric", "numeric", "numeric", 
                                      "numeric", "numeric"))

finalatmcsv   <-  read_excel("~/Documents/NR339/Climate/excel/finalatmcsv.xlsx", 
     col_types = c("numeric", "date", "numeric", 
         "numeric"))

rcndvi   <-  read.csv('/Users/simran/Documents/NR339/rcndvi.csv')

note: check your file paths and make sure to set your working directory to where you have stored your atmospheric and climate data

climate1  <-  data.frame(climateRC)
rcndvi  <-  data.frame(rcndvi)
atmdata  <-  data.frame(finalatmcsv)

Datasets

Let’s take look at the first few rows of the data frames:

Climate

kable(head(climateRC[1:5,]), align = "llcccccccr", format = "simple",  col.names = gsub("[.]", " ", names(climateRC)), digits = 3, format.args= list(scientifc = TRUE),
      caption = "Rock Creek Climate Data")
Rock Creek Climate Data
# Date Time Solar Radiation 3 Solar Radiation 4 Wind Speed, mph Gust Speed Temp, RH Rain Snow Depth Current
1 2019-09-04 17:30:00 114.4 43.1 0 2.26 81.509 28.8 NA -8.59 0.006
2 2019-09-04 18:00:00 148.1 30.6 0 2.26 81.554 29.1 0 -9.09 0.006
3 2019-09-04 18:30:00 101.9 20.6 0 2.82 79.788 31.3 0 -10.33 0.006
4 2019-09-04 19:00:00 34.4 6.9 0 0.00 75.160 39.6 0 -10.33 0.006
5 2019-09-04 19:30:00 0.6 0.6 0 0.00 71.017 48.9 0 -10.21 0.006

NDVI

kable(head(rcndvi[1:5,]), format = "simple",  align = "lcclr", digits = 3, format.args= list(scientifc = TRUE), caption = "NDVI Range Values")
NDVI Range Values
Date Count Minimum Maximum Range Mean
09-22-19 39 0.617 0.730 (0.6167,0.7297) 0.676
10-08-19 39 0.576 0.725 (0.5759,0.7251) 0.678
10-24-19 39 0.547 0.765 (0.5472,0.7653) 0.694
11-09-19 39 0.615 0.799 (0.615,0.7986) 0.701
11-25-19 39 0.650 0.758 (0.65,0.7581) 0.700

Atmospheric

kable(head(finalatmcsv[1:5,]), format = "simple", align = "lccr", col.names = c('Sample','Date', 'Absolute Pressure', 'Temperature (c)'), digits = 3, format.args = list (scientific = FALSE), caption = "Rock Creek Atmospheric Data")
Rock Creek Atmospheric Data
Sample Date Absolute Pressure Temperature (c)
1 2019-09-04 17:30:00 12.314 83.221
2 2019-09-04 18:00:00 12.313 86.823
3 2019-09-04 18:30:00 12.314 85.374
4 2019-09-04 19:00:00 12.321 78.274
5 2019-09-04 19:30:00 12.314 73.234

Looping

Now that we have these three data sets loaded in, we can begin to run through them to find our required data. Because NDVI data is recorded every 16 days, we will have to extract values from the other data frames over 16 day intervals.

Climate Data

If we look at our climate temperature data for the whole time series, it’s quite noisy:

alltemp<-plot(climate1$Date.Time, 
     climate1$Temp., main = "RC Temperature",
     xlab = "Date",  ylab = "Daily Temperature")

We have to clean this up a bit to get the values we need. So lets use this loop to run through all the climate data and extract mean temperature values for every 16 days in the data frame:

temphead  <-  head(climate1$Date.Time, n = 1)
temptail  <-  tail(climate1$Date.Time, n = 1)
starttemp = temphead
objecttemp = NULL
datet = NULL
while (starttemp <= temptail){
  endtemp= starttemp + ddays(16) 
  im  <-  data.frame(subset(climate1,
                        Date.Time >= starttemp & 
                          Date.Time <= endtemp))
  team   <-  mean(im$Temp)
  tempk = team 
  starttemp = endtemp 
  #taking each temp and storing it in object y which was null:
  objecttemp  <-  rbind(objecttemp,tempk)
  datet   <-  rbind(datet, starttemp)}

# date conversions for safety precautions
dates2  <-  as.Date.POSIXct(datet,
                         origin = "2019-09-04 17:30:00" )

We need to do the same thing for the min and max temperature:

dfhead  <-  head(climate1$Date.Time, n = 1) 
dftail  <-  tail(climate1$Date.Time, n = 1 )
start = dfhead
y   <-  NULL
x   <-  NULL

while (start <= dftail){
  end1 = start + ddays(16) 
  maxmin  <-  data.frame(subset(climate1, 
                              Date.Time >= start &
                                Date.Time <= end1))
  maxtempp  <-  max(maxmin$Temp.)
  mintempp  <-  min(maxmin$Temp.)
  mintemp = mintempp 
  maxtemp = maxtempp
  start = end1
#storing loop values of min/max temp in null object placeholder using rbind():
  y   <-  rbind(y, maxtemp)
  x   <-  rbind(x, mintemp)
}

Lets make a data frame of mean, min, and max temperature. Dont forget the dates!

temp  <-  data.frame("Mean Temp" = c(objecttemp), "Min Temp" = c(x), 
                    "Max Temp" = c(y), "Date" = dates2)
temp
Rock Creek Temperature (c)
Mean Temp Min Temp Max Temp Date
50.595 27.799 85.352 2019-09-20
45.156 19.929 81.730 2019-10-06
39.787 18.012 73.342 2019-10-22
37.275 8.577 70.286 2019-11-07
36.539 14.886 68.099 2019-11-23
28.858 4.813 52.029 2019-12-09

Now lets do some graphing!

 mg  <-  ggplot(data = temp, aes(x = dates2,y = objecttemp)) +
  geom_line(color = "lightblue")+ 
  geom_point(color = "lightblue")+
  xlab("dates") + ylab("temperature (C)") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + 
  scale_x_continuous(" ", labels = as.character(dates2),  breaks = dates2) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
        
finalmeant  <- mg + labs(
  title = " Rock Meadow Mean Temperature (C)", 
  subtitle = "June 2019 - July 2020")

finalmeant

Temperature decrease during winter months, and then rise to its highest in July of 2020.

Atmospheric Data

Now that we have our 20 climate data values for every 16 days within the time series, we should get our atmospheric pressure values in a similar format.

First we need to clean up and re-order our read in atmospheric data:

#removing temperature column from atmdata to make it easier to compile:
atmdata %>%   
  select(-c(temp))

#reformatting date to correlate with other climate data dates
atm   <-  atmdata[order(as.Date(atmdata$dated, format="%Y/%m/%d")),]

Now, using a similar loop to the climate data, we will run through our atmospheric data and extract the mean atmospheric pressure for every 16 days, and convert it from PSI to kPa:

#using loop to read thru atmospheric data to get mean pressure for 16 days:
atmhead  <-  head(atm$dated, n = 1)
atmtail  <-  tail(atm$dated, n = 1)
startatm = atmhead
atdf = NULL
dates1 = NULL

while (startatm <= atmtail){
  endatm= startatm + ddays(16) 
  atmdf  <-  data.frame(subset(atmdata, 
                           dated >= startatm & 
                             dated <= endatm))
  startatm = endatm 
  meanpres   <-  mean(atmdf$abspres) #taking mean pressure:
  kpa = meanpres * (6.89476)        #converting from psi to kpa
  atdf  <-  rbind(atdf, kpa)         #rbind() loop values to null object!
  dates1   <-  rbind(dates1, startatm)
}
#for some reason I had to reformat dates again and then make data frame 
atmosphere <- as.vector(atdf)
datesatm  <-  as.Date.POSIXct(dates1, origin = "2019-09-04 17:30:00" )
atmospheric   <-  data.frame("mean pressure (kpa)" = atmosphere, "dates" = datesatm)

Lets check out the data!

 at  <-  ggplot(data = atmospheric, aes(x = datesatm,y = atmosphere)) +
  geom_line()+ 
  xlab("Date") + ylab("Atmospheric Pressure (kPa)") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + 
  scale_x_continuous("Date", labels = as.character(datesatm), breaks = datesatm)
    
finalmeanatm  <- at + labs(
  title = " Rock Meadow Mean Atmospheric Pressure (kPa)", 
  subtitle = "June 2019 - July 2020") + 
  coord_cartesian(ylim = c(min(atdf), max(atdf)))

finalmeanatm

Looks pretty good! The range of atmospheric pressure is pretty small from 84.5 kPa - 85.25 kPa

NDVI

Now we will read in our NDVI data, which will be used in conjunction with the previously read in data to determine evapotranspiration. This data is already in a 16-day interval format, meaning no looping through is required. However we need to extract the min and max values, as well as the range.

ndvimax   <-  rcndvi[,"Maximum"]
ndvimin   <-  rcndvi[,"Minimum"]
# I converted it to numeric value and unlisted it using this:
ndvimax   <-  as.numeric(unlist(ndvimax))
ndvimin   <-  as.numeric(unlist(ndvimin))

# creating data frame of NDVI data and phi reference values
ndvidf1  <-  data.frame("NDVI Min"= ndvimin, "NDVI Max" = ndvimax, "α Min" = 0, 
                   "α Max" = 1.26, "Min Temp" = x, 
                   "Max Temp" = y, "Date" = rcndvi$Date)

Before I explain the interpolation scheme using NDVI, lets take a look at our NDVI data frame, and see what it entails so we can begin the next step of our calculation process.

Starting Values for α and NDVI min and max
NDVI Min NDVI Max α Min α Max Min Temp Max Temp Date
mintemp 0.617 0.730 0 1.26 27.799 85.352 09-22-19
mintemp.1 0.576 0.725 0 1.26 19.929 81.730 10-08-19
mintemp.2 0.547 0.765 0 1.26 18.012 73.342 10-24-19
mintemp.3 0.615 0.799 0 1.26 8.577 70.286 11-09-19
mintemp.4 0.650 0.758 0 1.26 14.886 68.099 11-25-19
mintemp.5 0.244 0.855 0 1.26 4.813 52.029 12-11-19
mintemp.6 0.251 0.575 0 1.26 7.279 48.094 12-27-19
mintemp.7 0.200 0.419 0 1.26 2.088 51.633 01-09-20
mintemp.8 0.275 0.657 0 1.26 4.379 44.173 01-25-20
mintemp.9 0.455 0.675 0 1.26 6.654 58.023 02-10-20
mintemp.10 0.482 0.706 0 1.26 13.771 61.032 02-26-20
mintemp.11 0.547 0.674 0 1.26 17.422 64.418 03-13-20
mintemp.12 0.420 0.565 0 1.26 2.012 57.722 03-29-20
mintemp.13 0.478 0.697 0 1.26 13.896 68.785 04-14-20
mintemp.14 0.500 0.635 0 1.26 23.518 79.743 04-30-20
mintemp.15 0.525 0.663 0 1.26 21.632 82.308 05-16-20
mintemp.16 0.543 0.702 0 1.26 25.201 87.939 06-01-20
mintemp.17 0.532 0.728 0 1.26 24.228 85.577 06-17-20
mintemp.18 0.535 0.735 0 1.26 31.892 95.760 07-03-20
mintemp.19 0.557 0.745 0 1.26 32.990 84.542 07-19-20

From my research these are the steps required to find phi/α (aerodynamic canopy resistance):

note that phi and α are interchangeable

  • Use Upper and lower bounds of α values to find α for each NDVI interval:

    • Global \(\ α_{min}\ \) = 0

    • Global \(\ α_{max}\ \) = 1.26

    • Phi min can be found for the Rock Creek NDVI interval using:

    • (RC NDVI min = 0, Global \(\ α_{min}\ \) = 0) & (RC NDVI max = max, Global \(α_{max}\) = 1.26)

    • α max is then calculated by using the lowest surface temperature pixel (\(T_{min}\)) within each interval:

      • (NDVI min = 0, α = RC \(\ α_{min}\ \) ) & (NDVI max = max, \(\ α_{max}\ \) = 1.26)
  • After bounds have been established for phi (α min, α max):
    • You linearly interpolate within each NDVI class between lowest temperature pixel:
    • (\(T_{min}\) , RC \(α_{max}\)) and ( \(T_{max}\), RC\(α_{min}\) )

Source: (PRiCE 1990), (Batra et al. 2006)

So…
Our first step is:

  • Use the two bounds where (RC NDVI min = 0, Global \(α_{min}\) = 0) & (RC NDVI max = max, Global \(α_{max}\) = 1.26) and the range of NDVI
    • (NDVI min, NDVI max) to interpolate α min for each NDVI interval.
    note: each step is explained within the code block in the code comments
# (NDVI MIN, α min), (NDVI = max, α MAX), and range
# two null values to store loop iteration values: 
apx   <-  NULL
daten  <-  NULL

#for every row in the ndvi data frame: 
for(row in 1:nrow(ndvidf1)) {
  
  # create objects g and f from each row and  column 1 and 2:
  g  <-  ndvidf1[row, 1]  # assigning object min ndvi 
  f  <-  ndvidf1[row,2]  #assigning object max ndvi :
  e  <-  ndvidf1[row,7]  #date from each row
  
  #concatenation time: 
  rangex  <-  c(0,g,f)            # range of x value is 0, tmin, tmax 
  rangey  <-  c(0,NA,1.26)        # range of y values is amin, NA, amax
                
  #the NA value is what we will be interpolating
  lp  <-  approx(rangey, rangex, xout = rangex, n = 3, method = "linear") 
  daten  <-  rbind(daten, e) #storing dates 
  lpy  <-  lp$y   
  final  <-  lpy[2]
  apx  <-  rbind(apx, final) #rbind() loop values to null object
}
ndvidf2  <-  data.frame("nmin"= ndvimin, "nmax" = ndvimax, "phi min" = apx,
                  "phi max" = 1.26, "mintemp" = x, 
                  "maxtemp" = y, "dated" = rcndvi$Date)

Now we have \(α_{min}\), which we can then use to interpolate α max for each NDVI interval:

# interpolation: (tmax, phimin) and (tmin, phimax) to find phimax
pm   <-  NULL
for(row in 1:nrow(ndvidf2)){
  # now we are using temperature to interpolate 
  tm  <-  ndvidf2[row, 5] # assigning min temp object 
  tx  <-  ndvidf2[row,6] # assigning max temp object 
  pmin  <-  ndvidf2[row,3]  # using interpolated phi min value 
  pmax  <-  1.26
  # temperature and phi values are now used to find phi max
  rangey  <-  c(tm, tx) # y values = min and max temp 
  rangex  <-  c(pmin, pmax) # x values = phi min and phi max 
  lp  <-  approx(rangex, rangey, n = 3, method = "linear") # linear interpolation 
  lpx  <-  lp$x
  final  <-  lpx[2]
  pm  <-  rbind(pm, final)}
ndvidf3 <-  data.frame("nmin"= ndvimin, 
                   "nmax" = ndvimax, 
                   "phi min" = apx, 
                   "phi max" = pm, 
                   "mintemp" =x, 
                   "maxtemp" = y, 
                   "dated" = rcndvi$Date)

Lastly, with our now interpolated \(α_{min}\) and \(α_{max}\), we can find the α of each NDVI interval:

 #find phi, which is between (tmax, phimin) and (tmin,phimx)
phi   <-  NULL
for(row in 1:nrow(ndvidf3))
  {
  tm  <-  ndvidf3[row, 5]
  tx  <-  ndvidf3[row,6]
  e  <-  ndvidf3[row,7]
  #using both interplated phi max and min values now to find phi 
  pmin  <-  ndvidf3[row,3]
  pmax  <-  ndvidf3[row,4]
  rangey  <-  c(tm, tx)
  rangex  <-  c(pmin, pmax)
  lp  <-  approx(rangex, rangey, n = 3, method = "linear") #linear interpolation
  #taking the 2nd x value of interpolation and establishing as phi
  lpx  <-  lp$x
  final  <-  lpx[2]
  #storing loop value into null placeholder object  using rbind()
  phi  <-  rbind(phi, final)}

In addition, Batra mentions their calculated interpolation scheme is represented in this equation:

\[α\ = \ 1.26⋅ \frac{T_{max} - T_{avg}}{∆_{temp}}\]

For the sake of comparing values for accuracy, we will now perform this calculation as well. The code is simple, and uses variables from the ndvidf data frame.

phicalc = 1.26*((y - objecttemp)/(y -x))
ndvidf4  <-  data.frame("Nmin"= ndvimin,"Nmax" = ndvimax, "α Min" = apx, 
                    "α Max" = pm, "α Interpolated" = phi,
                    "α Equation" = phicalc, "Min temp" = x, "Max Temp" = y,
                    "Date" = datecollected  <-  rcndvi$Date)
NDVI calculation values
Nmin Nmax α Min α Max α Interpolated α Equation Min temp Max Temp Date
final 0.617 0.730 0.357 0.809 0.583 0.761 27.799 85.352 09-22-19
final.1 0.576 0.725 0.331 0.796 0.564 0.746 19.929 81.730 10-08-19
final.2 0.547 0.765 0.332 0.796 0.564 0.764 18.012 73.342 10-24-19
final.3 0.615 0.799 0.390 0.825 0.607 0.674 8.577 70.286 11-09-19
final.4 0.650 0.758 0.391 0.826 0.608 0.747 14.886 68.099 11-25-19
final.5 0.244 0.855 0.165 0.713 0.439 0.618 4.813 52.029 12-11-19
final.6 0.251 0.575 0.115 0.687 0.401 0.566 7.279 48.094 12-27-19
final.7 0.200 0.419 0.066 0.663 0.365 0.600 2.088 51.633 01-09-20
final.8 0.275 0.657 0.143 0.702 0.422 0.440 4.379 44.173 01-25-20
final.9 0.455 0.675 0.244 0.752 0.498 0.657 6.654 58.023 02-10-20
final.10 0.482 0.706 0.270 0.765 0.518 0.729 13.771 61.032 02-26-20
final.11 0.547 0.674 0.292 0.776 0.534 0.767 17.422 64.418 03-13-20
final.12 0.420 0.565 0.188 0.724 0.456 0.640 2.012 57.722 03-29-20
final.13 0.478 0.697 0.264 0.762 0.513 0.714 13.896 68.785 04-14-20
final.14 0.500 0.635 0.252 0.756 0.504 0.741 23.518 79.743 04-30-20
final.15 0.525 0.663 0.277 0.768 0.523 0.729 21.632 82.308 05-16-20
final.16 0.543 0.702 0.303 0.781 0.542 0.714 25.201 87.939 06-01-20
final.17 0.532 0.728 0.308 0.784 0.546 0.695 24.228 85.577 06-17-20
final.18 0.535 0.735 0.312 0.786 0.549 0.671 31.892 95.760 07-03-20
final.19 0.557 0.745 0.330 0.795 0.562 0.721 32.990 84.542 07-19-20

Let’s get a visual representation of our α values, to see how they vary and how they change over our time series.

 plot(dates2, phi,type="o", col="blue", pch="o", lty=1, ylim= (0:1.26), 
      main = "α Interpolated vs α Equation", 
      xlab = "Date", ylab ="α")

# adding additional phi lines and creating legend for our plot 
points(dates2, phicalc, col="red", pch="*")
lines(dates2, phicalc, col="red",lty=2)
legend("topright",inset = 0.02, legend=c("α EQ", "α Interpolated"), 
       col=c("red", "blue"), lty=1:2, cex=0.8)

The α = 1.26equation values are slightly higher than the interpolated values, which makes sense as it purely used the global \(α_{max}\) instead of each NDVI interval’s \(α_{max}\) to determine α for each interval.

To compare accuracy and differences of the α values, lets use a paired t-test:

phistat  <-  t.test(
    phicalc,phi,alternative = c("two.sided", "less", "greater"), mu = 0,
    paired = TRUE, var.equal = FALSE, conf.level = 0.95)
phistat
## 
##  Paired t-test
## 
## data:  phicalc and phi
## t = 13.985, df = 19, p-value = 1.874e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1444800 0.1953361
## sample estimates:
## mean of the differences 
##                0.169908

With a p value of 1.874e-11 which is less than 0.5, and a t-statistic of 13.985, there is significant difference in mean values between the two α calculations. This means our α values are not alike, which means the using the global α as a constant creates higher α values of around 0.169908 or more for Rock Creek Meadow.

Psychometric Constant

Now that we have our climate, atmospheric, and NDVI calculated, we can determine the psychometric constant using atmospheric data.

The psychometric constant (γ) is equal to:

\[ y \ = \ C_p \ ⋅\ P ⋅ (\frac{R_a}{R_v}) \ ⋅ \ λ_m\]

\[ λ_m = \ Latent \ Heat \ of \ Vaporization \ = \ 2.435\frac{MJ}{Kg}\]

\[ C_p \ = \ Specific \ Heat \ of \ air \ (\frac{MJ}{\frac{kg}{C}})\]

\[ Where \ X \ =\ Humidity \ Ratio\ = \frac{R_a}{R_v} \ ⋅ \frac{E_a}{(\ P\ - \ E_a)}\]

\[ Where\ E_a\ = \ RH\ ⋅ \ E_s\]

\[Where\ E_s=\ Saturation \ vapor \ pressure \ of \ water \]

\[e_s = \ Saturation \ vapor\ pressure of water = (kPa)\]

\[e_a = \ Actual \ vapor\ pressure \ (kPa)\]

\[ P \ = \ Air \ Pressure\ (kPa)\]

\[\frac{R_a}{R_v}\ = \ Humidity \ Ratio\ (\frac{kg \ water\ vapor}{\ kg \ dry \ air})\ = \ 0.622\]

\[T_a = \ avg\ \ temperature \ (celsius)\]

There is a lot going on with the psychometric constant, so each section calculated will have its equation above for clarification.

Saturated Vapor Pressure, \(E_s\):

This is a combo of our previously determined variable:

\[ E_s \ = [1.007 +\ (3.46 \ ⋅ 10^{-5}\ ⋅ \ P)] ⋅ \ 6.1121 ⋅\ exp(\frac{(17.502 \ ⋅ \ T_a)}{(240.97 \ + \ T_a)}\]

#using temperature values and pressure values to find es 
esta   <- as.vector(objecttemp)
espres   <-  as.vector(atmosphere)

es = (1.007 + (3.46 * 10^-5 *espres))* 6.1121 * exp((17.502 * esta)/(240.97 * esta))

Relative Humidity and \(E_a\):

\[\ E_a\ = \ RH\ ⋅ \ E_s\]

We need to use a 16 day loop to find the average relative humidity values over 16 day intervals:

relativehumid   <-  climate1$RH
dfhead1  <-  head(climate1$Date.Time, n = 1)
dftail1  <-  tail(climate1$Date.Time, n = 1 )
started = dfhead
rh = NULL 
datesrh = NULL

while (started <= dftail1){
  end0 = started + ddays(16) 
  new1  <-  data.frame(subset
                   (climate1, 
                     Date.Time >= started & 
                     Date.Time <= end0))
  meanrh   <-  mean(new1$RH)  #taking mean relative humidity: 
  started = end0
  rh   <-  rbind(rh, meanrh) #storing loop values using rbind()
  datesrh   <-  rbind(datesrh, started)}
# Finding ES:
relativeh  <-  as.vector(rh)
ea = relativeh * es

Finding Humidity Ratio, X

\[\frac{R_a}{R_v}\ = \ Humidity \ Ratio\ (\frac{kg \ water\ vapor}{\ kg \ dry \ air})\ = \ 0.622\]

  • The ratio of \(R_a\) to \(R_v\) is ≈ 0.622:
    • \(R_a\):The gas constant for dry air= 286.9 \(\frac{J}{K}\)
    • \(R_v\): The gas constant for water vapor = 461.5 \(\frac{J}{K}\)
ex =(286.9/461.5) * ea/(espres - ea)

Finding Specific Heat of Air, \(C_p\)

\[C_p \ = \ Specific \ Heat \ of \ air \ (\frac{MJ}{\frac{kg}{C}}\ ) \ = \frac{1,005 + \ 1.82\ ⋅ X}{1,000}\]

cp = (1.005 + 1.82 * ex) / 1000
lm = 2.453

Finding the Psychometric Constant

\[y \ = \ C_p \ ⋅\ P ⋅ (\frac{R_a}{R_v}) \ ⋅ \ λ_m\]

psy = (cp * espres)/ ((286.9/461.5) * lm)
psycht <- data.frame('Psychometric Constant' = c(psy), 'Date' = c(dates2))
Psychometric Constant Date
-0.021 2019-09-20
-0.020 2019-10-06
-0.021 2019-10-22
-0.027 2019-11-07
-0.020 2019-11-23
-0.017 2019-12-09

Let’s take a look at our psychometric constant over our time series:

p   <-  ggplot(data = psycht, aes(x = dates2,y = psy)) +
  geom_line(color = "pink")+ 
  geom_point(color = "pink")+
  xlab("Dates") + ylab("Psychometric Constant") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + 
  scale_x_continuous("  ", labels = as.character(dates2), breaks = dates2)

finalp  <- p + labs(
  title = " Rock Meadow Mean Atmospheric Pressure (kPa)", 
  subtitle = "June 2019 - July 2020") + 
  coord_cartesian(ylim = c(-0.02679753, -0.01673537))
finalp

There is a drop in pressure on 11/07/19, which could be due to a drop in temperature, or data collection errors. If we line up temperature, pressure, and the psychometric constant, they all all drop significantly between 11/07/19 and 11/23/19.

Remember, we already graphed atmospheric pressure and temperature, so all we need to do is combine these three graphs into one area to visually compare:

at   <-  ggplot(data = atmospheric, aes(x = datesatm,y = atmosphere)) +
  geom_line(color = "pink")+ 
  geom_point(color = "pink")+
  xlab("Dates") + ylab("Atmospheric Pressure (kpa)") +
  theme(panel.grid.major = element_blank(),  panel.grid.minor = element_blank())
    
mg   <-  ggplot(data = temp, aes(x = dates2,y = objecttemp)) +
  geom_line(color = "lightblue")+ 
   geom_point(color = "lightblue")+
  xlab("Date") + ylab("Average Temperature (C)") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
 
p   <-  ggplot(data = psycht, aes(x = dates2,y = psy)) +
  geom_line(color = "orange")+ 
   geom_point(color = "orange")+
  xlab("Date") + ylab("Psychometric Constant") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
    
plot_grid(mg, at, p,
          ncol = 3, nrow = 1)

As displayed, there seems to be a simultaneous drop in pressure and the psychometric constant accompanied with a gradual drop in temperature. The temperature drop is more gradual due to a greater magnitude of change represented in temperature values compared to atmospheric pressure or the psychometric constant.

Now that we have the psychometric constant, we can begin to calculate the next variable in the Priestly-Taylor equation, which is solar radiation.

Solar Radiation, \(R_n\)

Solar Radiation values were provided in the climate csv in two columns. This is because the sensor that records solar radiation records two readings, solar radiation above, and solar radiation below. For the purpose of our calculations, these two values needed to be summed to find total solar radiation.

First, renaming the column names in our climate1 data frame will make it easier to manipulate:

climate1   <-  climate1 %>% rename(solarrad1 = Solar.Radiation...3)
climate1   <-  climate1 %>% rename(solarrad2 = Solar.Radiation...4)

Next, we will the same 16 day a loop iteration structure to obtain our solar radiation values over a 16 day interval within the time series. The \(R_n\) values were provided in \(\frac{\frac{Watts}{M^{2}}}{30 min}\) and need to be converted to \(\frac{\frac{MJ}{M^{2}}}{day}\)

dfhead  <-  head(climate1$Date.Time, n = 1)
dftail  <-  tail(climate1$Date.Time, n = 1)
start2 = dfhead
ss1   <-  NULL

while (start2 <= dftail){
  end2 = start2 + ddays(16) 
  new2  <-  data.frame(subset(climate1, 
                          Date.Time >= start2 &
                            Date.Time <= end2))
  solar1  <-  mean(new2$solarrad1)
  solar2  <-  mean(new2$solarrad2)
  start2 = end2
  solare  <-  solar1 + solar2  #take the sum of each interval, up + down = total 
  solared=solare *0.001
  solarcal=solared*3.6#converting from watts to MJ 
  ss1   <-  rbind(ss1, solarcal) #storing loop values using rbind():

}

Lets take a look at our solar radiation values over our time series:

solargraph  <-  data.frame("Solar Radiation" = ss1, 
                       "dates" = dates2)

solarg  <-  ggplot(data = solargraph, aes(x = dates2,y = ss1)) +
  geom_line(color = "darkblue")+ 
  geom_point(color = "darkblue")+
  xlab("Dates") + ylab("Average Solar Radiation ($$\frac{MJ}{m{^2}}/16 days$$") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
   theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
  scale_x_continuous("  ",labels = as.character(dates2), breaks = dates2)
solarg

This graph seems to follow our relative trend of lower values in the winter and highest values in July 2020, which correlates with the highest temperature values.

Now that we have our solar radiation values, we can find soil heat flux (G)

Finding Soil Heat Flux (G)

Soil heat flux is energy received by the soil to heat per unit of surface and time. If the soil is warmed, then soil heat flux is positive, if the soil is cooled, then the soil heat flux is negative

Soil heat flux can be calculated using a soil heat flux plate, or through remote sensing techniques, In the case of our data set, we did not have a soil heat flux plate at the study site, so we are using recommended equations using NDVI data and solar radiation to calculate soil heat flux. The equation is provided by Batra is as follows;

\[G \ = \ 0.583\ {exp \ (NDVI_a\ * R_n})\]

The NDVI value in this equation comes from the mean NDVI values for each interval which were provided in with the NDVI data, and \(R_n\) is solar radiation.

Let’s now calculate soil heat flux:

avgndvi   <-  rcndvi$Mean
sr   <-  as.vector(ss1)
g = 0.583 * exp(avgndvi * -2.31) * sr

Let’s take a look at g over our time series:

gdata  <-  data.frame("Soil Heat Flux" = g, 
                   "Dates" = dates2)

soilheat  <-  ggplot(data = gdata, aes(x = dates2,y = g)) +
  geom_line(color = "brown")+ 
  geom_point(color = "brown")+
  ylab("Soil Heat Flux") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5))+
  scale_x_continuous(" ",  labels = as.character(dates2),breaks = dates2)

soilheat

The data seems relatively consistent with our other variables, except there is an intense spike in late January and February of 2020, which we can perhaps attribute to the inconsistency of results due to this equation, our NDVI values, or solar radiation errors.

Calculating ∆

Saturated vapor pressure temperature curve :

\[∆ = \frac{4098*E_S}{T_a^{2}}\]

delta = 4098 * es / ((esta)^2)

Calculating Evapotranspiration

Now we have all our variables for calculation of actual evapotranspiration using the Priestly Taylor Equation!

Once again, the equation is:

\[ET = α \frac{( ∆ \ (Rn-G))} {(λ_v\ (\ ∆+γ \ ))}\]

\[α = \ Aerodynamic \ canopy \ resistance\] \[∆ = \ Slope \ of \ saturated \ vapor \ pressure \ temperature \ curve \ (\frac{kPa}{c})\] \[λv = \ Volumemetric \ latent\ heat\ of \ vaporization (\frac{2453\ MJ}{m^{3}})\] \[Rn =\ Net \ Radiation\ (\frac{MJ}{m^{2}}{day^{-1}})\] \[G = \ Soil \ Heat \ Flux\ Density \ (\frac{MJ}{m^{2}}{day^{-1}})\] \[γ = \ Psychometric \ constant (\frac{kPa}{C})\]

First find ET using the interpolated aerodynamic canopy resistance:

et1 = phi * delta *(ss1 - g)/ (lm *(delta + psy))

Now find ET with the provided α = 1.26 equation:

et2 = phicalc * delta *(ss1 - g)/(lm *(delta + psy))

#making data frame of the two et values 
evap  <-  data.frame("Equation ET" = c(et2), "Interpolated ET" = c(et1), 
                     "Date" = c(dates2))
Priestley-Taylor ET values for Rock Creek Meadow
Equation ET Interpolated ET Date
0.183 0.140 2019-09-20
0.155 0.117 2019-10-06
0.142 0.105 2019-10-22
0.091 0.082 2019-11-07
0.075 0.061 2019-11-23
0.035 0.025 2019-12-09
0.030 0.021 2019-12-25
0.047 0.028 2020-01-10
0.032 0.031 2020-01-26
0.084 0.064 2020-02-11
0.127 0.090 2020-02-27
0.134 0.093 2020-03-14
0.134 0.096 2020-03-30
0.216 0.155 2020-04-15
0.207 0.141 2020-05-01
0.197 0.141 2020-05-17
0.257 0.195 2020-06-02
0.277 0.218 2020-06-18
0.324 0.265 2020-07-04
0.239 0.186 2020-07-20

Let’s look at each ET individually over our time interval:

interpet  <-  ggplot(data = evap, aes(x = dates2,y = et1)) +
  geom_line(color = "blue")+ 
  geom_point(color = "blue")+
  xlab("Date") + ylab("Interpolated ET") +
  ggtitle("Rock Creek Meadow Interpolated Alpha ET")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + 
  scale_x_continuous("  ", labels = as.character(dates2), breaks = dates2)

calcet  <-  ggplot(data = evap, aes(x = dates2,y = et2)) +
  geom_line(color = "red")+ 
  geom_point(color = "red")+
  ggtitle("Rock Creek Meadow 1.26 Alpha ET")+
  xlab("Date") + ylab("1.26 alpha ET") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + 
  scale_x_continuous("  ", labels = as.character(dates2), breaks = dates2)

The Interpolated α ET looks like:

interpet

The 1.26 α ET looks like:

calcet

By just visual analysis, these follow a similar trend of ET increasing in the summer months. In addition, they seem to be around the same range of values. This means our two different methodologies for calculating ET provide similar results.

Using Evapotranspiration package to find Prietsly Taylor ET

Now, using one more method to find Priestly Taylor ET will be performed using the R Evapotranspiration Package.

Similar loops will be used, to find daily values for ET, which we will then use and take the mean value of 16 day to compare to our calculated values.

To learn more about the package, check out these links:

R Evapotranspiration Documentation
CRAN Package ‘Evapotranspiration’

There are various functions in this package which use different methodologies to determine ET, for our case we will be using the ET.PriestleyTaylor function.

The variables required for this function include:

  • \(T_{max}\) = Temperature Max (celsius)
  • \(T_{min}\) = Temperature Min (celsius)
  • \(RH\) = Relative Humidity
  • \(RH_{min}\) = Relative Humidity Min
  • \(RH_{max}\) = Relative Humidity Max
  • \(R_s\) = Solar Radiation
  • \(V_s\) = Saturated Vapor Pressure
  • \(U_z\) = Wind speed (mph)
  • Precip = rain (inches)
  • Year
  • Month
  • Days

These are the correct titles for the objects we read into the ET function, and they are name and case specific to ensure proper calculations.

The only thing that has changed in our loop is that now we are storing all our values, instead of 16 days, we must do 1 day intervals, as the package only performs sub daily, daily or monthly calculations. Once we have our daily ET value, we will find the mean value for every 16 days.

First, we need to load all our loops in again, so lets do that. By now, this looping process should make sense. I also re read our climate data for this, just to make sure that it was the raw data.

climateRC   <-  read_excel("~/Documents/NR339/Climate/excel/climateRC.xls",
                           col_types = c("numeric", "date", "numeric",
                                         "numeric", "numeric", "numeric",
                                         "numeric", "numeric", "numeric",
                                         "numeric", "numeric"))
## New names:
## * `Solar Radiation` -> `Solar Radiation...3`
## * `Solar Radiation` -> `Solar Radiation...4`
atmrc <- read_excel("~/Documents/NR339/Climate/excel/finalatmcsv.xlsx",
                    col_types = c("skip", "date", "numeric","skip"))

#------------------------Temp, Tmin, Tmax-------------------------------------

#Tmin and Tmax 
dfhead  <-  head(climateRC$`Date Time`, n = 1)
dftail  <-  tail(climateRC$`Date Time`, n = 1)
start = dfhead
a   <-  NULL
b  <-  NULL
while (start <= dftail){ 
  end1 = start +ddays(1)
  maxmin  <-  data.frame(subset(climateRC, 
                                `Date Time` >= start &
                                  `Date Time` <= end1))
  mintemp =  min(maxmin$Temp.) 
  maxtemp = max(maxmin$Temp.)
  start = end1
  a   <-  rbind(a, maxtemp)
  b  <-  rbind(b, mintemp)
}


Tmin<-as.vector(b)
Tmax<-as.vector(a)

#Temp - mean temp per day 
temphead<-head(climateRC$`Date Time`, n = 1)
temptail<- tail(climateRC$`Date Time`, n = 1)
starttemp = temphead
tt = NULL
while (starttemp <= temptail){
  endtemp= starttemp +ddays(1)
  im<-data.frame(subset(climateRC,
                        `Date Time` >= starttemp & 
                          `Date Time` <= endtemp))
  team <- mean(im$Temp.) 
  starttemp = endtemp 
  tt<- rbind(tt,team)
  }

Temp<-as.vector(tt)

#----------------------------------Rh, Rhmin, Rhmax-----------------------------

#Used one loop to calculate all three: 
dfheadRH<-head(climateRC$`Date Time`, n = 1)
dftailRH<- tail(climateRC$`Date Time`, n = 1)
startedRH = dfheadRH
rh = NULL 
rhmi = NULL
rhm = NULL
while (startedRH <= dftailRH){
  endRH = startedRH + ddays(1)
  newRH<-data.frame(subset(climateRC, `Date Time` >= startedRH 
                           & `Date Time` <= endRH))
  startedRH = endRH
  rhday<-mean(newRH$RH)
  rhmaxim<-max(newRH$RH)
  rhminim<-min(newRH$RH)
  rh <- rbind(rh, rhday)
  rhmi <- rbind(rhmi, rhminim)
  rhm<- rbind(rhm, rhmaxim)
}

RHmax<-as.vector(rhm)
RHmin<-as.vector(rhmi)
RH<-as.vector(rh)

#-----------------------------Solar Radiation,Rs--------------------------------

dfheadss<-head(climateRC$`Date Time`, n = 1)
dftailss<- tail(climateRC$`Date Time`, n = 1)
startsolar = dfheadss
ss1 <- NULL
dates<-NULL
while (startsolar <= dftailss){
  end2 = startsolar + ddays(1)
  new2<-data.frame(subset(climateRC, `Date Time` >= startsolar & `Date Time` <= end2))
  s1<-sum(new2$Solar.Radiation...3)
  s2<-sum(new2$Solar.Radiation...4)
  dat<-as.Date(end2)

  convert<- (s1 + s2) * 0.001
  Rs <- convert / 3.6 #converting from watts to MJ
  startsolar = end2
  dates<-rbind(dates, dat)
  ss1 <- rbind(ss1, Rs)
}

#storing climate dates for future 'Days', 'Months', 'Year'
date1<-as.Date(dates, origin = "1970-01-01, 0:00:00 ")
Rs <- ss1

#------------------------------------Precip-------------------------------------

#rain 
dfhead1<-head(climateRC$`Date Time`, n = 1)
dftail1<- tail(climateRC$`Date Time`, n = 1)
start2 = dfhead
inch <- NULL
while (start2 <= dftail){
  end2 = start2 + ddays(1)
  new2<-data.frame(subset(climateRC, `Date Time` >= start2
                          & `Date Time` <= end2))
  water<-sum(new2$`Rain`, 
             na.rm = TRUE) 
  #there is an Na value in our dataset, we dont want that
  mm= water *25.4    #converting from inches to MM 
  start2 = end2
  inch <- rbind(inch, mm)
}

Precip<-as.vector(inch)

#-------------------------------Wind speed, Uz---------------------------------- 

dfhead1<-head(climateRC$`Date Time`, n = 1)
dftail1<- tail(climateRC$`Date Time`, n = 1)
start2 = dfhead
mph <- NULL
while (start2 <= dftail){
  end2 = start2 + ddays(1)
  new2<-data.frame(subset(climateRC, `Date Time` >= start2 
                          & `Date Time` <= end2))
  wind<-mean(new2$Wind.Speed..mph)
  start2 = end2
  mph <- rbind(wind, mph)
}
uz <- as.vector(mph)

#--------------------------Saturated Vapor Pressure, V--------------------------

dfhead1<-head(atmrc$dated, n = 1)
dftail1<- tail(atmrc$dated, n = 1)
startatm = dfhead
psi = NULL
while (startatm <= dftail1){
  endatm= startatm + ddays(1)
  atmdf<-data.frame(subset(atmrc, dated >= startatm & dated <= endatm))
  startatm = endatm 
  pres <- mean(atmdf$abspres)
  kpa = pres * (6.89476) #converting from psi to kpa
  atdf<- rbind(psi, kpa)
}

kpapres <- as.vector(atdf)

#using the same equation in earlier equations to find Vs
Vs = (1.007 + (3.46 * 10^-5 *kpapres))* 6.1121 * exp((17.502 * Temp)/(240.97 * Temp))

#------------------------------------Formatting date---------------------------------
Year = lubridate::year(date1)
Month = lubridate::month(date1)
Day = lubridate::day(date1)
Hour = lubridate::hour(date1)

There are also constants required for the calculation, which must be stored in list format. These constants in include:
* Elevation (ft)

  • Lambda (latent heat of vaporization) = \((\frac{2.453\ MJ}{m^{3}})\)

  • Latitude in radians (lat_rad)

  • Solar constant (Gsc) =\((\frac{0.0820 MJ}{m^{2}}{min^{-1}})\)

  • Stefan-Boltzman constant (Sigma) \(4.903*10^{-9} \frac{MJ}{K^{4}}/m^{2}/day\)

  • Soil heatflux (G negligible for daily time step) = 0

  • Aerodynamic canopy resistance (which is α or AlphaPt) = 1.26

\[ET = α \frac{( ∆ \ (Rn-G))} {(λ_v\ (\ ∆+γ \ ))}\]

\[α = \ Aerodynamic \ canopy \ resistance\] \[∆ = \ Slope \ of \ saturated \ vapor \ pressure \ temperature \ curve \ (\frac{kPa}{c})\]

\[λ_v = \ Volumemetric \ latent\ heat\ of \ vaporization (\frac{2.453\ MJ}{m^{3}})\] \[R_n =\ Net \ Radiation\ (\frac{MJ}{m^{2}}{day^{-1}})\]

\[G = \ Soil \ Heat \ Flux\ Density \ (\frac{MJ}{m^{2}}{day^{-1}})\] \[γ = \ Psychometric \ constant (\frac{kPa}{C})\]

constants <- list(Elev = 5100, 
                  lambda = 2.45, 
                  lat_rad =  0.3547, 
                  Gsc = 0.0820, 
                  sigma = (4.903*10^-9), 
                  alphaPT = 1.26,
                  G = 0)

Now that we have determined all our variables and constants, now use the ReadInputs() function from the Evapotranspiration package to prep our data for calculation. This requires creating a data frame with your variables and then using it in the ReadInputs function. Click here for an in depth description of the ReadInputs() function. For our purposes, we want these the function parameters to be:

  • varnames = c("Tmax", "Tmin", "RHmin", "RHMax". "Vs", "uz", "Precip", "Rs", "Year", "Month","Day"

  • stopmissing = c(1,1,1)

  • timestep = "daily"

  • interp_missing_entries = TRUE

  • missing_method = "neighbouring_average"

  • abnormal_method= "neighbouring_average"

  • interp_abnormal = TRUE

  • message = "yes"

processdata<-data.frame(Tmax, Tmin, Temp, RHmin, RHmax,RH, Vs, uz, Precip, Rs, Year, Month,Day,Hour)

preist <- ReadInputs(varnames = c("Tmax", "Tmin","Vs", "uz", "Precip", 
          "RHmin", "RHmax", "Rs","Vs", "Year","Month","Day"),
          processdata, constants, 
          stopmissing = c(1,1,1), 
          timestep = "daily", interp_missing_days = FALSE, 
          interp_missing_entries = FALSE,
          missing_method = "neighbouring average", interp_abnormal = TRUE,
          abnormal_method = "neighboring average",
          message = "yes")
## The maximum acceptable percentage of date indices is 1 %
## The maximum acceptable percentage of missing data is 1 %
## The maximum acceptable percentage of continuous missing data is 1 %

Now once we have read our climate and atmospheric inputs, they can be processed using the function ET.PriestlyTaylor for final calculation’s.

Click here here for the documentation for the ET.PriestlyTaylor function. We will be using our recently interpreted preist input data, the constants list and:

  • ts (time setting) = "daily"

  • solar (\(R_s\)) = "data"

  • alpha (Surface albedo for Rock Creek) = 0.65

  • message = "yes"

  • AdditionalStats = "yes"

finalpt <- ET.PriestleyTaylor(preist, 
                        constants, 
                        ts = "daily",  
                        solar = "data", 
                        alpha = 0.65, 
                        message = "yes", 
                        AdditionalStats =  "yes")
## Priestley-Taylor Potential ET
## Evaporative surface: user-defined, albedo = 0.65
## Solar radiation data have been used directly for calculating evapotranspiration
## Timestep: daily
## Units: mm
## Time duration: 2019-09-05 to 2020-07-06
## 306 ET estimates obtained
## Basic stats
## Mean: 0.48
## Max: 1.17
## Min: -1.06
#-----------------------Obtaining results for ET over 16 day interval-----------

results<-finalpt$ET.Daily

dailyAET<-fortify.zoo(results)
dailyAET 
pthead<-head(dailyAET$Index, n = 1)
pttail<- tail(dailyAET$Index, n = 1)
ptstart = pthead
ptet = list()
ptdates = list()

while(ptstart <= pttail){
  ptend= ptstart + ddays(16)
  endtime<-data.frame(subset(dailyAET,
                             Index >= ptstart & 
                               Index <= ptend))
  etp<-mean(endtime$results)
  dated = ptend
  ptet[[length(ptet)+1]]= etp
  ptdates[[length(ptdates)+1]] = dated
  ptstart=ptend
}

unlisted<-unlist(ptdates)
finalpriestly<-data.frame("Priestly ET" = unlist(ptet), "Date" = as.Date(unlisted))

Now lets take a look at the finalpriestly results:

kable(finalpriestly, format = "simple", col.names=gsub("[.]", " ", names(finalpriestly)), align = "lcr",
     digits = 3, format.args = list(scientific = FALSE), caption = "Package Priestley-Taylor ET values for Rock Creek Meadow")

The values are slightly higher compared two our other calculation methods, we can graph the results to see if they follow a similar increasing trend.

plotetp  <-  ggplot(data = finalpriestly, aes(x = Date,y = Priestly.ET)) +
  geom_line(color = "green")+ 
  geom_point(color = "green")+
  xlab("Date") + ylab(" ET Package Priestly AET ") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + 
  scale_x_continuous("  ", labels = as.character(dates2), breaks = dates2)
    
finalplotted  <- plotetp + labs(
  title = "Rock Meadow R Package Evapotranspiration", 
  subtitle = "June 2019 - July 2020") + 
  coord_cartesian(ylim = c(min(unlist(ptet)), max(unlist(ptet))))

finalplotted

With this data, there is more variation between high and low values, and there seems to be a general trend increasing values until July 20th, where it drops significantly. However, the overall trend of the increase in ET in the summer months follows our other data sets as well. Lets compare all three graphs now:

interpet1  <-  ggplot(data = evap, aes(x = dates2,y = et1)) +
  geom_line(color = "blue")+ 
  geom_point(color = "blue")+
  xlab("Date") + ylab("Interpolated ET") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

calcet1  <-  ggplot(data = evap, aes(x = dates2,y = et2)) +
  geom_line(color = "red")+ 
  geom_point(color = "red")+
  xlab("Date") + ylab("1.26 alpha ET") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

packageet  <- ggplot(data = finalpriestly, aes(x = Date,y = Priestly.ET)) +
  geom_line(color = "green")+ 
  geom_point(color = "green")+
  xlab("Date") + ylab(" ET  Package") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) 
    
plot_grid(calcet1, interpet1, packageet, align = "hv", axis = "l",
         ncol = 3, nrow = 1, lables = c("1.26 alpha ET", "interpolated ET", "R Package ET"), label_size = 10, label_x = "Date", 
         label_y = "ET" )
## Warning in as_grob.default(plot): Cannot convert object of class character into
## a grob.

This follows a similar trend of increasing evapotranspiration in the summer months, then with a steep drop in late July. However, it seems that the ET.PriestlyTaylor() method’s values are larger overall, as the Y-axis ranges from 0 to 0.8. Lets overlap all of these values to get a sense of how similar or not similar they are.

 plot(dates2, ptet ,type="o", col="green", pch="o", lty=1, ylim=c(0,1), 
      main = "Priestly-Taylor AET of Rock Creek ", 
      xlab = "Date", ylab ="ET (mm/16 days)")

# adding additional phi lines and creating legend for our plot 
points(dates2, et1, col="blue", pch="*")
lines(dates2, et1, col="blue",lty=2)
points(dates2, et2, col ="red", pch="-")
lines(dates2, et2, col = "red", lty=2)
legend("topleft",inset = 0.02, legend=c("Package ET", "ET Calculated", "ET Interpolated"), 
       col=c("green","red", "blue"), lty=1:2, cex=0.8)

As previously mentioned, the Evapotranspiration package ET values are at a much larger scale than our values that were calculated manually. This might represent discrepancies in the calculation method and data collections. In this case, it serves us best to compare and contrast our manually calculated values since they are at the same scale, and then perhaps try another method of the Evapotranspiration method to compare our modeled values to.

We can use one more calculation method from the ET package to find ET of Rock Creek Method. Specifically, we can run the standard Penman-Montieth calculation model to find another set of ET values to compare our package values to. For simplification, the Pennman-Montieth Method is very similar to the Priestly Taylor method, but the Priestly Taylor method has simplified the vapor deficit and convection terms are reduced to a single empirical constant (α). The equation is listed as:

\[λ_vET_{p,om}= \frac{ Δ⋅\ (R_n−G)\ +ρ_a⋅c_p⋅(\frac{e_s−e_a}{r_a})}{Δ + γ⋅(1 + \frac {r_s}{r_a})}\]

\[ET_p = Penman \ Monteith \ Evapotranspiration \ (\frac{m}{day})\]

\[λ_v = \ Volumetric\ latent \ heat \ of \ vaporiazation\ (\frac{2543 \ MJ}{m^3})\]

\[Δ = \ Slope \ of \ the \ saturation\ vapor \ pressure-temperature \ curve\ (\frac{kPa}{C})\]

\[R_n = \ Net \ radiation \ (\frac{\frac{MJ}{m^2}}{day})\]

\[G = \ Soil\ heat \ flux \ density\ (\frac{\frac{MJ}{m^2}}{day})\]

\[ρ_a = \ Air\ density \ for \ a \ given \ air \ pressure \ (\frac{kg}{m^3})\]
\[c_p = \ Specific \ heat\ of\ air \ (\frac{\frac{MJ}{kg}}{C})\]

\[e_s = \ Saturation \ vapor \ pressure (kPa)\]

\[e_a = \ Actual \ vapor\ pressure \ (kPa)\]

\[r_a = \ Air \ resistance \ (\frac{day}{m})\]

\[ γ = \ Psychometric \ constant (\frac{kPa}{C})\]

\[r_s = \ Surface \ resistance (\frac{day}{m})\]

Now lets use the Evapotranspiration package’s ET.Penman() function to determine ET of Rock Creek Meadow once more. We are using ET.Penman() because our data fits the function requirements which include:

Tmax, Tmin, RHmax, RHmin, Rs

processdata<-data.frame(Tmax, Tmin, Temp, RHmin, RHmax,RH, Vs, uz, Precip, Rs, Year, Month,Day,Hour)

constants <- list(Elev = 5000, 
                  lambda = 2.45, 
                  lat_rad =  0.3547, 
                  Gsc = 0.0820, 
                  sigma = (4.903*10^-9), 
                  G = 0,
                  z = 2) 

penman1 <- ReadInputs(varnames = c("Tmax", "Tmin","Vs", "uz", "Precip", 
          "RHmin", "RHmax", "Rs","Vs", "Year","Month","Day"),
          processdata, constants, 
          stopmissing = c(1,1,1), 
          timestep = "daily", interp_missing_days = FALSE, 
          interp_missing_entries = FALSE,
          missing_method = "neighbouring average", interp_abnormal = TRUE,
          abnormal_method = "neighboring average",
          message = "yes")
## The maximum acceptable percentage of date indices is 1 %
## The maximum acceptable percentage of missing data is 1 %
## The maximum acceptable percentage of continuous missing data is 1 %
penman2<-ET.Penman(penman1, constants, 
                            ts = "daily",  
                            solar = "data",
                            alpha = 0.65, 
                            wind = "yes",
                            message = "yes", 
                            windfunction_ver = 1948, 
                            AdditionalStats =  "yes")
## Penman Potential ET
## Evaporative surface: user-defined, albedo = 0.65 ; roughness height = 0.001 m
## Solar radiation data have been used directly for calculating evapotranspiration
## Wind data have been used for calculating the Penman evaporation. Penman 1948 wind function has been used.
## Timestep: daily
## Units: mm
## Time duration: 2019-09-05 to 2020-07-06
## Basic stats
## Mean: 1.69
## Max: 4.14
## Min: 0.07
#-----------------------Obtaining results for ET over 16 day interval-----------
penmanresult<-penman2$ET.Daily
penmanAET<-fortify.zoo(penmanresult)

thead<-head(penmanAET$Index, n = 1)
ttail<- tail(penmanAET$Index, n = 1)
tstart = thead
tet = list()
tdates = list()

while(tstart <= ttail){
  tend= tstart + ddays(16)
  pendtime<-data.frame(subset(penmanAET,
                             Index >= tstart & 
                               Index <= tend))
  pen<-mean(pendtime$penmanresult)
  pdated = tend
  tet[[length(tet)+1]]= pen
  tdates[[length(tdates)+1]] = pdated
  tstart=tend
}
unlistedpen<-unlist(tdates)
finalpenman<-data.frame("Penman ET" = unlist(tet), "Date" = as.Date(unlistedpen))

Lets look at our ETPenman() results:

Package Penman ET values for Rock Creek Meadow
Penman ET Date
2.023 2019-09-21
1.608 2019-10-07
2.197 2019-10-23
2.650 2019-11-08
1.821 2019-11-24
0.671 2019-12-10
0.574 2019-12-26
1.161 2020-01-11
0.586 2020-01-27
1.587 2020-02-12
2.004 2020-02-28
1.698 2020-03-15
1.148 2020-03-31
1.875 2020-04-16
1.862 2020-05-02
1.857 2020-05-18
1.894 2020-06-03
2.187 2020-06-19
2.659 2020-07-05
2.412 2020-07-21

Now using the ETComparison() function we can compare and contrast ET.PriestlyTaylor() and ETPenman():

ETComparison(penman2, finalpt, type = "Daily", ylim=c(-1,8),
labs=c("RC Penman"," RC Priestly Taylor"))

The non-exceedance probability is the probability of an uncertain parameter exceeding a certain threshold. ETPenman() has a higher exceedance probability. ETPenman() overall range and values are slightly higher than ET.PriestlyTaylor().

Lets now graph the two 16 day interval Penman and Priestly Taylor values to compare:

 plot(dates2, tet ,type="o", col="purple", pch="o", lty=1, ylim=c(0,5), 
      main = "Average Penman vs Priestly Taylor ET of Rock Creek ", xlab = "Date", ylab ="ET (mm/day)")

# adding additional phi lines and creating legend for our plot 
points(dates2, ptet, col="green", pch="*")
lines(dates2, ptet, col="green",lty=2)
legend("topleft",inset = 0.02, legend=c("Penman ET", "Prieslty ET"),
       col=c("purple","green"), lty=1:2, cex=0.8)

The Priestly Taylor values differ significantly from the Penman values, but they seem to still increase over the summer months and drop in the winter.

Okay now that we have all our data, lets compare and do some exploring!

Analysis

interpet1  <-  ggplot(data = evap, aes(x = dates2,y = et1), color = "blue") +
  geom_line(color ="blue")+ 
  geom_point(color ="blue")+
  xlab("Date") + ylab("Interpolated ET") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

calcet1  <-  ggplot(data = evap, aes(x = dates2,y = et2), color = "red") +
  geom_line(color = "red")+
  geom_point(color = "red")+
  xlab("Date") + ylab("1.26 alpha ET") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 


packageet  <- ggplot(data = finalpriestly, aes(x = Date,y = Priestly.ET), color = "green") +
  geom_line(color = "green")+ 
  geom_point(color = "green")+
  xlab("Date") + ylab("Package ET Priestly") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) 


penmanet  <- ggplot(data = finalpenman, aes(x = Date,y = Penman.ET), color = "purple") +
  geom_line(color = "purple")+ 
  geom_point(color = "purple")+
  xlab("Date") + ylab("Package ET Penmman ") +
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank()) 

plot_grid(packageet,calcet1, penmanet, interpet1)

Based off of the data, overall overall evapotranspiration of Rock Creek Meadow decreased from November 2019 to January of 2020, and then began increasing as temperatures increased up until July 2020. To test how much of an influence temperature was on our evapotranspiration data, lets explore the data.

finalone<-data.frame("Penman ET" = unlist(tet), "Package Priestly ET" = unlist(ptet), "ET Interpolated" = et1, "ET Alpha 1.26" = et2, "Date" = dates2, "Temp" = objecttemp, "Alpha" = phi)
finalone
Compiled Calculated ET Values for Rock Creek Meadow
Penman ET Package Priestly ET ET Interpolated ET Alpha 1 26 Date Temp Alpha
final 2.023 0.189 0.140 0.183 2019-09-20 50.595 0.583
final.1 1.608 0.150 0.117 0.155 2019-10-06 45.156 0.564
final.2 2.197 0.540 0.105 0.142 2019-10-22 39.787 0.564
final.3 2.650 0.692 0.082 0.091 2019-11-07 37.275 0.607
final.4 1.821 0.371 0.061 0.075 2019-11-23 36.539 0.608
final.5 0.671 0.416 0.025 0.035 2019-12-09 28.858 0.439
final.6 0.574 0.398 0.021 0.030 2019-12-25 29.761 0.401
final.7 1.161 0.568 0.028 0.047 2020-01-10 28.025 0.365
final.8 0.586 0.408 0.031 0.032 2020-01-26 30.274 0.422
final.9 1.587 0.652 0.064 0.084 2020-02-11 31.256 0.498
final.10 2.004 0.785 0.090 0.127 2020-02-27 33.687 0.518
final.11 1.698 0.639 0.093 0.134 2020-03-14 35.794 0.534
final.12 1.148 0.624 0.096 0.134 2020-03-30 29.403 0.456
final.13 1.875 0.736 0.155 0.216 2020-04-15 37.671 0.513
final.14 1.862 0.394 0.141 0.207 2020-05-01 46.669 0.504
final.15 1.857 0.452 0.141 0.197 2020-05-17 47.180 0.523
final.16 1.894 0.291 0.195 0.257 2020-06-02 52.363 0.542
final.17 2.187 0.521 0.218 0.277 2020-06-18 51.723 0.546
final.18 2.659 0.434 0.265 0.324 2020-07-04 61.752 0.549
final.19 2.412 0.030 0.186 0.239 2020-07-20 55.058 0.562

Lets look at the density of each ET methodology:

etden1<- (density(finalone$ET.Interpolated, bw = "nrd0"))

d1<-plot(etden1,
     main  = "ET Interp Density")

etpen1<- (density(finalone$Penman.ET, bw = "nrd0"))

d2<-plot(etpen1,
     main  = "ET Penman Density")

etalphad<- (density(finalone$ET.Alpha.1.26, bw = "nrd0",))

d3<-plot(etalphad, 
     main  = "ET Alpha 1.26 Density")

etpackpt<- (density(finalone$Package.Priestly.ET, bw = "nrd0"))
d4<- plot(etpackpt,
     main  = "ET Package Priestly Density")

There seems to be an uneven distribution of values for the Penman Fortieth and the Priestly Taylor Evapotranspiration Method. The NDVI methodology seems to have similar densities.

te<-ggplot(data = finalone, aes(x = Date))+
  xlab("Date") + ylab("ET (mm / 16 day)") +
  ggtitle("Rock Creek ET and Avg Temperature")+
  geom_point(mapping = aes(y = et1, color = objecttemp), color = "blue")+
  geom_smooth(mapping = aes(y = et1, color = objecttemp), color = "blue", method = "loess")+
  geom_point(mapping = aes(y = et2, color = objecttemp), color = "red")+
  geom_smooth(mapping = aes(y = et2, color = objecttemp), color = "red", method = "loess")+
  geom_point(mapping = aes(y = Penman.ET, color = objecttemp), color = "green")+
  geom_smooth(mapping = aes(y = Penman.ET, color = objecttemp), color = "green", method = "loess")+
  geom_point(mapping = aes(y = Package.Priestly.ET, color = objecttemp), color = "purple")+
  geom_smooth(mapping = aes(y = Package.Priestly.ET, color = objecttemp), color = "purple", method = "loess")
  
tempvset<- te + labs(colour = "Temperature (C)", show.legend = TRUE)

at<-ggplot(data = finalone, aes(x = Date))+
  xlab("Date") + ylab("Phi") +
  ggtitle("Rock Creek ET and Alpha")+
  geom_point(mapping = aes(y = phi, color = et1), color = "blue")+
  geom_smooth(mapping = aes(y = et1, color = et1), color = "blue")+
  geom_point(mapping = aes(y = phi, color = et2), color = "red")+
  geom_smooth(mapping = aes(y = et2, color = et2), color = "red")+
  geom_point(mapping = aes(y = phi, color = Penman.ET), color = "green")+
  geom_smooth(mapping = aes(y = phi, color =Penman.ET), color = "green")+
 geom_point(mapping = aes(y = Package.Priestly.ET, color = Package.Priestly.ET), color= "purple")+
  geom_smooth(mapping = aes(y = phi, color = Package.Priestly.ET), color="purple")

alphat<- at + labs(color = "ET", show.legend = TRUE)

alphat 
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

tempvset
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

Here we have overlaid each ET method with temperature or α values to show which values of correspond with which calculation method.
The lower ET values correlate to lower temperature values, where as the higher ET values correlate to higher temperature values. This is also seen in our α values, where there seems to greater distribution from the mean with higher α values with higher temperatures and ET values.

Conclusion

Overall, there seems to be a trend in higher temperature and α values leading to higher levels of ET, which makes sense scientifically due to higher temperatures causing greater evaporation.

When determining which method was most effective for calculating ET, it seems that the using NDVI data provided similar results, regardless of using a constant α value or not. Using the Evapotranspiration package, it seems that the results are in a closer range to each other when compared to the NDVI data.

One issue with the methodology might be that because Batra’s method was meant for large scale areas, whereas Rock Creek Meadow is a much smaller area than what was used in the original study. This is why NDVI values might make ET values smaller for a smaller area.

Further research is required to determine how successful NDVI data is at determining evapotranspiration of a small area.

References

Evapotranspiration documentation

Batra, Namrata, Shafiqul Islam, Virginia Venturini, Gautam Bisht, and LE Jiang. 2006. “Estimation and Comparison of Evapotranspiration from Modis and Avhrr Sensors for Clear Sky Days over the Southern Great Plains.” Remote Sensing of Environment 103 (1): 1–15.

Jiang, Le, and Shafiqul Islam. 2001. “Estimation of Surface Evaporation Map over Southern Great Plains Using Remote Sensing Data.” Water Resources Research 37 (2): 329–40.

PRiCE, John C. 1990. “Using Spatial Context in Satellite Data to Infer Regional Scale Evapotranspiration.” IEEE Transactions on Geoscience and Remote Sensing 28 (5): 940–48.