This document is written to describe the calculation procedure of crop evapotranspiration using the dual crop coefficient approach. The data used in this article is from High-Tunnel study performed by the author and the others1 in 2018. There is source code containing global parameters and functions to execute some of the codes below. These resources will be released as a seperate R-script file in the future when the package development for this research project is accomplished. The script file contains many useful calculators, constants, unit converters, and user created functions. Most of the FAO-56 equations (Allen et al. 1998) needed for the dual crop coefficient method are incorporated into this file. To learn more about basic concepts and principles on crop evapotranspiration, consult the FAO-56.

Introductinon

There are a number of ways to calculate or estimate crop evapotranspiration (ETc) (Colaizzi et al. 2017; S. R. Evett et al. 2012; Shaughnessy et al. 2015). Among these the FAO-56’s Penman-Monteith method (Allen et al. 1998) is based on measurements of meteorological variables that affect ETc in combination with measurements of crop and soil characteristics, i.e. crop coefficients, height of crop, and daily soil water balance (see below).

The single crop coefficient (Kc) approach begins with a calculation of reference crop evapotranspiration (ET0) from the measured meteorological factors followed by a calculation of Kc over the growing season. Kc can be referred to known values reported in the FAO-56 (Table 12 for Kc) and one can use the referred values to construct Kc curves. The initial Kc curves then can be further adjusted by climatic factors and height of crop. The last step is to multiply ET0 by Kc.

\(ET_c = ET_0 * K_c\)

The dual Kc approach uses the same ET0 value, but breaks Kc further down into two sub-elements – the basal crop coefficient (Kcb) and the transpiration coefficient (Ke).

\(K_c = K_{cb} + K_e\)

Kcb can be referred to the known values in the same FAO-56 (Table 18 for Kcb). Likewise, Kcb curves over the growing season should be constructed and adjusted like the same way used in the Kc. Ke calculation relys on daily calculation of soil water balance. Recording of water input to soil from precipitaion and irrigation, therefore, becomes an integral part in this calculation. All wetting events should be tracked and provided to estimate Ke and this also depends on the characteristics of the soil such as total evaporable water (TEW), readily evaporable water (RAW), field capacity (FC), wilting point (WP), etc. The type or composition of the soil heavily affects these values. The calculation of daily Ke is technical, requiring a computational tool. All information about these metrics are given and fully detailed in the FAO-56. The abbridged methodology paper is in composition as of March 21, 2019 ( OAP_ET_WUE_Cal, to be updated).

\(ET_c = ET_0 * (K_{cb} + K_e)\)

The first component, ET0 Kcb_ explains transpirational water loss on the crop surfaces while the second component, ET0 Ke explains evaporational water loss on the soil surfaces.

Procedure

The entire procedure can be summarized in the following steps.

  1. Import packages & source code
  2. Import & tidy data sets – weather data (met), phenology data (time), irrigation data (irg), & yield data (yd)
  3. Calculate ET0 using met
  4. Import & adjust Kc (Kcb) using met & time
  5. Construct Kc (Kcb) curves
  6. Calculate Ke & ETc adding irg to the available data sets created
  7. Verify the results
  8. Calculate WUE using yd & the available data sets created

Step 1: Import packages & source codes

library(tidyverse)  # includes many useful data manipulation & exploration pkgs
library(lubridate)  # deals with timeseries data
library(DT)  # creates interactive data tables
library(scales)  # enables graphical modifications on ggplot2 objects
library(directlabels)  # enables graphical modifications on ggplot2 objects
library(nlme)  # builds & analyzes linear mixed effect model 
library(lsmeans)  # conducts a post hoc comparison 
library(knitr)  # compiles html objects & generates a report file
library(openxlsx)  # exports tabulated dataframes or R-objects to xlsx files

source("0_ht_src.R")  # provides useful unit conversion & calculation tools

The source code also contains package loading functions that are outlined above.

Step 2: Import & tidy data sets

The meteorological data acquisition was performed by two sets of weather stations. Set 1 was a pair of two weather stations installed on site. For Set 1, a compact data logger ( CR300, Campbell Sci. Inc., Logan, UT was hooked up with various weather sensors to obtain climate variables. Set 2 was a weather station ( RX3000, Onset Computer Co., Bourne, MA) standing in close proximity of the experimental plots. The data from Sets 1 and 2 were combined, reorganized, and compiled into a single data set for convenient further calculation processes. Due to a large size of the data, the raw data is not presented here, but rather the processed, summarized data will be presented later in Step 3.

The phenology & physiology data were collected by various instruments and tools. Only height of the crops will be used in the further steps. I extract hgt after summarizing the raw data, time.

The irrigation data were corrected manually by Jimmy Gray. The irrigation records will be used in the calcuation of daily soil water balance for Ke cal.

Step 3: Calculate ET0

The FAO-56 is implemented into the CalET0 function. Details can be found in the source code.

The following is to back up some missing data points in the weather data set.

et0 is set aside to a seperate data frame for future convenience.

Step 4: Import & adjust Kc (Kcb)

Then, the Kcb table from the FAO-56 comes. The lengths of each stage (denoted as dat in kcbcrv below) were determined by the observed phenology data above. A total of five K-values is needed – “st”, “ini”, “mid”, “end”, and “tm”.

Prepare some materials by slicing the stages of the growth & development.

Since doy is set up now, calculations for average h can follow. Additionally, average ws and min_rh are required to adjust Kcb values.

There are missing two averages of hgt in HT data sets. These are omitted during joining and tidying the weather data sets due to several outliers. The outliers were resulted from the power outage problem in the tunnels. We need to fix these before going further.

Using the averages of h, ws, and min_rh, complete the adjustments to Kcb values.

Step 5: Construct Kc (Kcb) curves

Let’s check out the adjusted Kcb values in a graph. This can be considered Kcb curves except the data points expressed as the lines on the graph represent theoretical values.

These graphs above are from just five points of Kcb values. Actual Kcb values should be produced to further calculate Ke values on a daily basis. To do this, functions to connect the adjusted kcvals during the crop development stage (from “ini” to “mid”) and the late-season stage (from “end” to “tm”).

Simple straight line models allow for extrapolating K-values during development and termination of the growth.

All the materials are ready to produce the Kcb curve data set.

Compile the calculated and tidied data into one data frame.

Let’s compare the full Kcb curves generated by the functions and the previously built Kcb curves using the five referenced points to see if there were computational errors.

Step 6. Calculate Ke & ETc

Step 6.1. Initializing

Bring all the data sets for the calculation together. We need et0 and kcbfull data frames.

Step 6.2. Iterating

Below is the function implemented in the iteration process to adjust and calculate correct Dest, Dend, Kr, Ke, E, and other parms.

for (i in seq_along(kekcb$stg)) {  # for loop cals
  if (kekcb$dest[i] <= bl_rew) {
      kekcb$stg[i] <- 1
      kekcb$kr[i] <- 1
      kekcb$ke[i] <- CalKe(kekcb$kr[i], kekcb$kmax[i], kekcb$kcb[i], kekcb$few[i])
      kekcb$e[i] <- kekcb$et0[i] * kekcb$ke[i]
      kekcb$kc[i] <- kekcb$ke[i] + kekcb$kcb[i]
      kekcb$etc[i] <- kekcb$et0[i] * kekcb$kc[i] 
      if (i == 1) {
        kekcb$dend[i] <- bl_tew + CalDei(kekcb$dest[i], kekcb$prec[i], kekcb$irwtr[i], kekcb$fw[i], kekcb$e[i], kekcb$few[i], 0, 0)
        kekcb$dest[i+1] <- kekcb$dend[i]
      } else {
        kekcb$dend[i] <- min(max(CalDei(kekcb$dest[i], kekcb$prec[i], kekcb$irwtr[i], kekcb$fw[i], kekcb$e[i], kekcb$few[i], 0, 0), 0), bl_tew)  # min & max setups to ensure 0 <= dend <= tew
        if (i < 91) {
          kekcb$dest[i+1] <- kekcb$dend[i] 
          }
      }
    } else {
      kekcb$stg[i] <- 2
      kekcb$kr[i] <- CalKr(bl_tew, bl_rew, kekcb$dest[i])
      kekcb$ke[i] <- CalKe(kekcb$kr[i], kekcb$kmax[i], kekcb$kcb[i], kekcb$few[i])
      kekcb$e[i] <- kekcb$et0[i] * kekcb$ke[i]
      kekcb$kc[i] <- kekcb$ke[i] + kekcb$kcb[i]
      kekcb$etc[i] <- kekcb$et0[i] * kekcb$kc[i]
      kekcb$dend[i] <- min(max(CalDei(kekcb$dest[i], kekcb$prec[i], kekcb$irwtr[i], kekcb$fw[i], kekcb$e[i], kekcb$few[i], 0, 0), 0), bl_tew)  # min & max setups to ensure 0 <= dend <= tew
      if (i < 91) {
        kekcb$dest[i+1] <- kekcb$dend[i] 
        }
    } 
}

Step 6.3. Compiling

Cumulative T, E, ETc, along with cumulative total water input (I+P) are calculated here too.

Step 7. Verify the results

Kc curves are constructed by the dual Kc method for each treatment.

Using this information, ETc graphs are constructed with separation of ETc into E (black shadeds) and T (green shades) components.

Better comparisons can be achieved by plotting ETc on the same plotting space for each irrigation management.

Cumulative water loss by ETc can be plotted in the same way.

Cumulative stats can be summarized in a table.

Step 8. Calculate WUE

This is an extra step to make calculations of water use efficiency of the crops in each treatment group. Water use efficiency can be defined by the following equation.

\(WUE = \frac{YD}{\sum ET_c}\)

The yield data is also required to further calculate water use efficiency of the crops with each treatment.

Statistical analysis is conducted on the raw data to test the effects of irrigation method (noted irg) on calculated WUE. A linear mixed-effect model was used to run two-way ANOVA for testing sys, spp, and the interaction effects. sys and spp were set to a fixed variable and zn (Zone) was set to a random variable in the model. The Tukey’s method was used in post hoc comparisons between treatments.

md.wue1 <- lme(data = wue, wue ~ sys * spp, random = ~1 | zn)
anova(md.wue1)  # no interaction effects, drop the irg:spp from the model
##             numDF denDF   F-value p-value
## (Intercept)     1     9 155.08176  <.0001
## sys             1     9  48.10337  0.0001
## spp             1     9  50.65120  0.0001
## sys:spp         1     9   0.31558  0.5880
md.wue2 <- lme(data = wue, wue ~ sys + spp, random = ~1 | zn)
anova(md.wue2)  # p < .001 for sys*** & .0021 for spp**
##             numDF denDF   F-value p-value
## (Intercept)     1    10 155.08177  <.0001
## sys             1    10  51.63753  <.0001
## spp             1    10  54.37255  <.0001
summary(md.wue2)  
## Linear mixed-effects model fit by REML
##  Data: wue 
##        AIC      BIC    logLik
##   45.12955 47.95429 -17.56477
## 
## Random effects:
##  Formula: ~1 | zn
##         (Intercept)  Residual
## StdDev:   0.4737926 0.6641837
## 
## Fixed effects: wue ~ sys + spp 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  6.020212 0.3726038 10 16.157140       0
## sysOF       -2.386387 0.3320918 10 -7.185926       0
## sppTom      -2.448770 0.3320918 10 -7.373775       0
##  Correlation: 
##        (Intr) sysOF 
## sysOF  -0.446       
## sppTom -0.446  0.000
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.7354646 -0.4606874  0.1248534  0.2376666  1.5782902 
## 
## Number of Observations: 16
## Number of Groups: 4
lsmeans(md.wue2, pairwise ~ sys, adjust = "tukey")  # ht:of = a:b 
## $lsmeans
##  sys lsmean    SE df lower.CL upper.CL
##  HT    4.80 0.334  3     3.73     5.86
##  OF    2.41 0.334  3     1.35     3.47
## 
## Results are averaged over the levels of: spp 
## d.f. method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE df t.ratio p.value
##  HT - OF      2.39 0.332 10 7.186   <.0001 
## 
## Results are averaged over the levels of: spp
lsmeans(md.wue2, pairwise ~ spp, adjust = "tukey")  # pep:tom = a:b
## $lsmeans
##  spp lsmean    SE df lower.CL upper.CL
##  Pep   4.83 0.334  3     3.77     5.89
##  Tom   2.38 0.334  3     1.32     3.44
## 
## Results are averaged over the levels of: sys 
## d.f. method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast  estimate    SE df t.ratio p.value
##  Pep - Tom     2.45 0.332 10 7.374   <.0001 
## 
## Results are averaged over the levels of: sys
pwue <- wue %>% 
  filter(spp == "Pep")

twue <- wue %>% 
  filter(spp == "Tom")

md.p1 <- lme(data = pwue, wue ~ sys, random = ~1 | zn)
anova(md.p1)  # ANOVA table, p = .0084 for sys**
##             numDF denDF  F-value p-value
## (Intercept)     1     3 156.8842  0.0011
## sys             1     3  21.3514  0.0191
summary(md.p1)  # regression model summary
## Linear mixed-effects model fit by REML
##  Data: pwue 
##      AIC      BIC    logLik
##   25.925 25.09204 -8.962502
## 
## Random effects:
##  Formula: ~1 | zn
##         (Intercept) Residual
## StdDev:   0.6072962 0.671212
## 
## Fixed effects: wue ~ sys 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  5.923567 0.4525854  3 13.088285  0.0010
## sysOF       -2.193097 0.4746185  3 -4.620758  0.0191
##  Correlation: 
##       (Intr)
## sysOF -0.524
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.1912211 -0.5234339  0.0621916  0.4200628  1.2315562 
## 
## Number of Observations: 8
## Number of Groups: 4
lsmeans(md.p1, pairwise ~ sys, adjust = "tukey")  # ht:of = a:b
## $lsmeans
##  sys lsmean    SE df lower.CL upper.CL
##  HT    5.92 0.453  3     4.48     7.36
##  OF    3.73 0.453  3     2.29     5.17
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE df t.ratio p.value
##  HT - OF      2.19 0.475  3 4.621   0.0191
md.t1 <- lme(data = twue, wue ~ sys, random = ~1 | zn)
anova(md.t1)  # ANOVA table, p = .0072 for sys**   
##             numDF denDF  F-value p-value
## (Intercept)     1     3 57.88957  0.0047
## sys             1     3 39.21393  0.0082
summary(md.t1)  # regression model summary
## Linear mixed-effects model fit by REML
##  Data: twue 
##        AIC      BIC   logLik
##   23.81908 22.98612 -7.90954
## 
## Random effects:
##  Formula: ~1 | zn
##         (Intercept)  Residual
## StdDev:   0.4702277 0.5825862
## 
## Fixed effects: wue ~ sys 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  3.668087 0.3743397  3  9.798819  0.0023
## sysOF       -2.579677 0.4119507  3 -6.262103  0.0082
##  Correlation: 
##       (Intr)
## sysOF -0.55 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -0.97387247 -0.42879952 -0.05708827  0.31546487  1.52240300 
## 
## Number of Observations: 8
## Number of Groups: 4
lsmeans(md.t1, pairwise ~ sys, adjust = "tukey")  # ht:of = a:b
## $lsmeans
##  sys lsmean    SE df lower.CL upper.CL
##  HT    3.67 0.374  3    2.477     4.86
##  OF    1.09 0.374  3   -0.103     2.28
## 
## d.f. method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE df t.ratio p.value
##  HT - OF      2.58 0.412  3 6.262   0.0082

To conclude the data analysis, let’s take a look at WUE of the crops under different irrigation management practices. The two-way ANOVA results and following post hoc analysis results are marked on the graphs. Statistical significance codes are provided at P < .10, .05, .01, .001 with o, *, **, and ***, respectively. Different alphabet letters indicate significant differences at P < .05 level.

References

Allen, Richard G., Luis S. Pereira, Dirk Raes, and Martin Smith. 1998. “Crop Evapotranspiration (guidelines for computing crop water requirements).” Rome: FAO. doi:10.1016/j.eja.2010.12.001.

Colaizzi, Paul D., Susan A. O’Shaughnessy, Steve R. Evett, and Ryan B. Mounce. 2017. “Crop evapotranspiration calculation using infrared thermometers aboard center pivots.” Agricultural Water Management 187. Elsevier B.V.: 173–89. doi:10.1016/j.agwat.2017.03.016.

Evett, Steven R., Robert C. Schwartz, Terry A. Howell, R. Louis Baumhardt, and Karen S. Copeland. 2012. “Can weighing lysimeter ET represent surrounding field ET well enough to test flux station measurements of daily and sub-daily ET?” Advances in Water Resources 50. Elsevier Ltd: 79–90. doi:10.1016/j.advwatres.2012.07.023.

Shaughnessy, Susan A O, Steven R. Evett, Paul D. Colaizzi, Susan A. O’Shaughnessy, Steven R. Evett, and Paul D. Colaizzi. 2015. “Dynamic prescription maps for site-specific variable rate irrigation of cotton.” Agricultural Water Management 159. Elsevier B.V.: 123–38. doi:10.1016/j.agwat.2015.06.001.

  1. Dr. Paul Colaizzi and Melanie Baxter from USDA-CPRL at Bushland, Drs. Charles Rush, Qingwu Xue from Texas A&M AgriLife Research at Amarillo, James Gray, Jewel Arthur, Jared Bull, student workers from Texas A&M AgriLife Research at Bushland.↩

 

A work by Tony Rho

hyungmin.rho@tamu.edu

Creative Commons License