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 the Ogallala Aquifer Program 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.

Introduction

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 13, 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

source("0_oap_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 were collected by a weather station ( RX3000, Onset Computer Co., Bourne, MA). 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.

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

Step 3: Calculate ET0

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

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

## Error: $ operator is invalid for atomic vectors

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.

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

## # A tibble: 9 x 15
## # Groups:   irg, spp [9]
##   irg   spp   grp   yd_avg  yd_sd  yd_n  yd_se wue_avg wue_sd wue_n wue_se
##   <chr> <chr> <chr>  <dbl>  <dbl> <int>  <dbl>   <dbl>  <dbl> <int>  <dbl>
## 1 ctrl  corn  ctrl~  1.03  0.147      3 0.0851   11.9   1.69      3  0.978
## 2 mch   corn  mch-~  0.981 0.251      3 0.145    11.9   3.04      3  1.76 
## 3 pvt   corn  pvt-~  0.641 0.0546     3 0.0315    6.78  0.578     3  0.334
## 4 ctrl  pep   ctrl~  2.64  0.137      3 0.0790   24.1   1.25      3  0.722
## 5 mch   pep   mch-~  3.64  0.605      3 0.349    36.0   5.98      3  3.45 
## 6 pvt   pep   pvt-~  1.12  0.245      3 0.142     9.43  2.06      3  1.19 
## 7 ctrl  tom   ctrl~  2.91  0.0739     3 0.0427   23.6   0.599     3  0.346
## 8 mch   tom   mch-~  3.75  2.39       3 1.38     32.5  20.7       3 12.0  
## 9 pvt   tom   pvt-~  1.76  0.635      3 0.367    13.1   4.73      3  2.73 
## # ... with 4 more variables: ieff_avg <dbl>, ieff_sd <dbl>, ieff_n <int>,
## #   ieff_se <dbl>

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 irg, spp, and the interaction effects. irg 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 ~ irg * 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    12 169.65914  <.0001
## irg             2     6  11.70594  0.0085
## spp             2    12   8.91672  0.0042
## irg:spp         4    12   1.60077  0.2373
md.wue2 <- lme(data = wue, wue ~ irg + spp, random = ~1 | zn)
anova(md.wue2)  # p = .0108 for irg* & .0038 for spp**
##             numDF denDF   F-value p-value
## (Intercept)     1    16 152.95217  <.0001
## irg             2     6  10.55321  0.0108
## spp             2    16   8.03866  0.0038
summary(md.wue2)  
## Linear mixed-effects model fit by REML
##  Data: wue 
##        AIC      BIC    logLik
##   177.2799 184.9172 -81.63997
## 
## Random effects:
##  Formula: ~1 | zn
##          (Intercept) Residual
## StdDev: 0.0004411243 7.903084
## 
## Fixed effects: wue ~ irg + spp 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  11.242505  3.400946 16  3.305700  0.0045
## irgmch        6.924005  3.725549  6  1.858519  0.1125
## irgpvt      -10.093698  3.725549  6 -2.709318  0.0351
## spppep       13.001570  3.725549 16  3.489840  0.0030
## spptom       12.871092  3.725549 16  3.454817  0.0033
##  Correlation: 
##        (Intr) irgmch irgpvt spppep
## irgmch -0.548                     
## irgpvt -0.548  0.500              
## spppep -0.548  0.000  0.000       
## spptom -0.548  0.000  0.000  0.500
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.53812650 -0.36240718 -0.02825509  0.34532627  2.69800496 
## 
## Number of Observations: 27
## Number of Groups: 9
lsmeans(md.wue2, pairwise ~ irg, adjust = "tukey")  # ctrl:mch:pvt = a, a, b
## $lsmeans
##  irg  lsmean   SE df lower.CL upper.CL
##  ctrl  19.87 2.63  8    13.79     25.9
##  mch   26.79 2.63  6    20.34     33.2
##  pvt    9.77 2.63  6     3.33     16.2
## 
## Results are averaged over the levels of: spp 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast   estimate   SE df t.ratio p.value
##  ctrl - mch    -6.92 3.73  6 -1.859  0.2304 
##  ctrl - pvt    10.09 3.73  6  2.709  0.0783 
##  mch - pvt     17.02 3.73  6  4.568  0.0091 
## 
## Results are averaged over the levels of: spp 
## P value adjustment: tukey method for comparing a family of 3 estimates
lsmeans(md.wue2, pairwise ~ spp, adjust = "tukey")  # corn:pep:tom = b, a, a
## $lsmeans
##  spp  lsmean   SE df lower.CL upper.CL
##  corn   10.2 2.63  6     3.74     16.6
##  pep    23.2 2.63  6    16.74     29.6
##  tom    23.1 2.63  6    16.61     29.5
## 
## Results are averaged over the levels of: irg 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast   estimate   SE df t.ratio p.value
##  corn - pep   -13.00 3.73 16 -3.490  0.0080 
##  corn - tom   -12.87 3.73 16 -3.455  0.0086 
##  pep - tom      0.13 3.73 16  0.035  0.9993 
## 
## Results are averaged over the levels of: irg 
## P value adjustment: tukey method for comparing a family of 3 estimates
cwue <- wue %>% 
  filter(spp == "corn")

pwue <- wue %>% 
  filter(spp == "pep")

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

md.c1 <- lme(data = cwue, wue ~ irg, random = ~1 | zn)
anova(md.c1)  # ANOVA table, p = .0337 for irg*
##             numDF denDF   F-value p-value
## (Intercept)     1     6 224.97679  <.0001
## irg             2     6   6.28693  0.0337
summary(md.c1)  # regression model summary
## Linear mixed-effects model fit by REML
##  Data: cwue 
##        AIC      BIC    logLik
##   38.86256 37.82136 -14.43128
## 
## Random effects:
##  Formula: ~1 | zn
##         (Intercept)  Residual
## StdDev:    1.907577 0.7153414
## 
## Fixed effects: wue ~ irg 
##                 Value Std.Error DF   t-value p-value
## (Intercept) 11.881883  1.176232  6 10.101650  0.0001
## irgmch       0.013614  1.663443  6  0.008184  0.9937
## irgpvt      -5.101439  1.663443  6 -3.066795  0.0220
##  Correlation: 
##        (Intr) irgmch
## irgmch -0.707       
## irgpvt -0.707  0.500
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -0.51724158 -0.06127470 -0.01336973  0.11489006  0.53061131 
## 
## Number of Observations: 9
## Number of Groups: 9
lsmeans(md.c1, pairwise ~ irg, adjust = "tukey")  # ctrl:mch:pvt = a, a, b
## $lsmeans
##  irg  lsmean   SE df lower.CL upper.CL
##  ctrl  11.88 1.18  8     9.17    14.59
##  mch   11.90 1.18  6     9.02    14.77
##  pvt    6.78 1.18  6     3.90     9.66
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast   estimate   SE df t.ratio p.value
##  ctrl - mch  -0.0136 1.66  6 -0.008  1.0000 
##  ctrl - pvt   5.1014 1.66  6  3.067  0.0501 
##  mch - pvt    5.1151 1.66  6  3.075  0.0496 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
md.p1 <- lme(data = pwue, wue ~ irg, random = ~1 | zn)
anova(md.p1)  # ANOVA table, p < .001 for irg***
##             numDF denDF  F-value p-value
## (Intercept)     1     6 349.7712  <.0001
## irg             2     6  38.3796   4e-04
summary(md.p1)  # regression model summary
## Linear mixed-effects model fit by REML
##  Data: pwue 
##        AIC      BIC    logLik
##   46.08614 45.04493 -18.04307
## 
## Random effects:
##  Formula: ~1 | zn
##         (Intercept) Residual
## StdDev:    3.482667    1.306
## 
## Fixed effects: wue ~ irg 
##                 Value Std.Error DF   t-value p-value
## (Intercept)  24.13896  2.147448  6 11.240763  0.0000
## irgmch       11.85100  3.036951  6  3.902269  0.0080
## irgpvt      -14.70534  3.036951  6 -4.842141  0.0029
##  Correlation: 
##        (Intr) irgmch
## irgmch -0.707       
## irgpvt -0.707  0.500
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -0.58138235 -0.11203728 -0.01115692  0.12319419  0.54500903 
## 
## Number of Observations: 9
## Number of Groups: 9
lsmeans(md.p1, pairwise ~ irg, adjust = "tukey")  # ctrl:mch:pvt = b, a, c
## $lsmeans
##  irg  lsmean   SE df lower.CL upper.CL
##  ctrl  24.14 2.15  8    19.19     29.1
##  mch   35.99 2.15  6    30.74     41.2
##  pvt    9.43 2.15  6     4.18     14.7
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast   estimate   SE df t.ratio p.value
##  ctrl - mch    -11.9 3.04  6 -3.902  0.0187 
##  ctrl - pvt     14.7 3.04  6  4.842  0.0069 
##  mch - pvt      26.6 3.04  6  8.744  0.0003 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
md.t1 <- lme(data = twue, wue ~ irg, random = ~1 | zn)
anova(md.t1)  # ANOVA table, ns, ctrl:mch:pvt = a, a, a
##             numDF denDF  F-value p-value
## (Intercept)     1     6 31.69909  0.0013
## irg             2     6  1.87063  0.2337
summary(md.t1)  # regression model summary
## Linear mixed-effects model fit by REML
##  Data: twue 
##        AIC      BIC    logLik
##   60.42437 59.38317 -25.21218
## 
## Random effects:
##  Formula: ~1 | zn
##         (Intercept) Residual
## StdDev:     11.5035 4.313811
## 
## Fixed effects: wue ~ irg 
##                  Value Std.Error DF   t-value p-value
## (Intercept)  23.579336  7.093175  6  3.324229  0.0159
## irgmch        8.907402 10.031265  6  0.887964  0.4087
## irgpvt      -10.474311 10.031265  6 -1.044167  0.3366
##  Correlation: 
##        (Intr) irgmch
## irgmch -0.707       
## irgpvt -0.707  0.500
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -0.61469802 -0.01959149  0.01208837  0.06081302  0.56797760 
## 
## Number of Observations: 9
## Number of Groups: 9

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. Drs. Paul Colaizzi, Qingwu Xue, & Charles Rush↩