Crop Evapotranspiration Calculation by Soil Water Balance Method
Tony Rho, PhD, Texas A&M AgriLife Research at Amarillo
Published on April 05, 2019
Updated on April 16, 2019
This document is written to describe the calculation procedure of crop evapotranspiration using the soil water balance 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 completed. The script file contains many useful calculators, constants, unit converters, and user created functions.
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 one considered most accurate is using a mass balance or a soil water balance method. Although both methods incorporate relatively heavy deployment of instruements compared with other commercially competitive methods, the accuracy excels since the loss of water by the plant and the soil surfaces is basically measured, not estimated like the method based on meteorological measurement of weather variables ( FAO-56)(Allen et al. 1998).
The soil water balance method requires measurements of soil water content (SWC) in any forms at given time intervals. The SWC measurement can be done with either manual data collection by a neutron probe or automatic data collection by an intelligent sensor. The measurement should be made at different layers of the soil by varying the depth of the measurement. With the collected soil water profile, one can estimate how much of water is consumed by crops and at which layer(s) the water is taken. The data collection process is often laborious and technically difficult. So, many researchers rely on automated measurement by field sensors. Typically, field sensors can be deployed directly to the soil at varying depths where crop(s) of interest stands. When combined with measurements of water fluxes to the designated soil volume (a.k.a. the control volume), water loss by the crop stand – ETc – can be calculated using the following equation that is used to calculate ETc.
\(ET_c = I + P + F + R + \Delta S\)
where I is irrigation, P is precipitation, F is net subsurface flux into the control volume, R is net runoff or runon to the control volume surface, and \(\Delta S\) is the net change in volumetric SWC in the control volume. All of the measures are in ‘mm’ units. In many cases of agricultural lands with flat soil surfaces, R is assumed to be 0.
Concurrently, the crop coefficient – Kc – can be estimated by solving the equation below.
\(K_c = \frac{ET_c}{ET_0}\)
where Kc is crop coefficient, ETc and ET0 stand for crop evapotranspiration and reference crop evapotranspiration each.
Procedures
- Import packages & source code
- Import & tidy data sets – weather data (
met), irrigation data (irg), soil water data (swc), & yield data (yd) - Calculate ET0 using
met - Calculate ETc using the SWC equation
- Calculate Kc from ETc
- Calculate WUE using
yd& the available data sets created
Steps 1-3 are a duplicate of those in Crop Evapotranspiration Calculation by Dual Crop Coefficient Method. The calculations presented in Steps 4 & 5 are unique to this report.
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 toolsThe 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.
The volumetric SWC data were corrected by a set of time-domain reflectometer (TDR) probes on each treatment plot. The TDR probes ( TDR-315L, Acclima Inc., Meridian, ID) were installed at .1, .2, .3, .45, .7, and 1 m depths of the control volume of the soil. SWC was measured at 5-min time intervals. Hourly SWC dataframe was based on these measurements. Daily SWC dataframe was subsampled from hourly SWC by extracting SWC at the beginning of the day (swcini, 00:00:00). \(\Delta S\) is calculated by subtracting SWC at the end of the day (swcend, 23:55:00) from SWC at the beginning of the day.
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. Calculate ETc using the SWC equation
Step 4.1. Calculate \(\Delta S\)
Let’s combine et0 and swc to compile the weather data and the soil water data, wtrbal.
The compiled dataframe contains et0, etc, swc, and kc to plot over time.
See the below table and find more details about the properties of the soil type. PAW/WP/FC stands for plant available water/wilting point/field capacity of the soil. The raw data measured by the TDR probes are presented without any modifications. Note that there are some values more than the field capacity (FC) of the soils. FC and WP of the Pullman Soils at ~ 1 m depth is 357 and 187 mm, respectively. These are indicated in select plots down below.
Skimming through hourly SWC can provide an insight as to how the soil responds to input or output of water.
Hourly SWC on Jul 20, 2018 shows a typical response to a wetting event.
Contrary to the Jul 20, hourly SWC on Jul 21 shows a typical response to a drying event without any water input.
Now, let’s hourly overview the development of SWC over the season.
SWC on a daily basis can be plotted along with wetting events – irrigation and precipitation.
To better compare the differences in SWC responses between HT and OF, the SWC are plotted on the days after transplanting (DAT) scale. There was a 3 week difference in the transplanting dates being the transplanting in HT three weeks earlier.
Harvesting was made 123 DAT for OF vegetables, meaning that the data for HT vegetables should be modified to have the same length. There is data points up to 146 DAT for HT vegetables. To parallel the differences in the estimates, the last 23 days of the data for HT should be dropped from the dataframe.
Step 4.2. Calculate ETc & Cumulatives
ETc can be calculated using day-to-day \(\Delta S\) of this graph. Following is the soil water balanced ETc and cumulative ETc over the season. Plots are on the DAT scale.
Cumulative ETc and water input from irrigation plus precipitation can be plotted separately.
Step 5. Calculate Kc from ETc
Kc is calculated by dividing ETc from the soil water balance method by ET0 from the meteorological method.
However, the time-course of Kcs seems incorrect for some reason. There are too many outliers to make a conclusion on the results. Kcs should fall into a 0-1.5 range. Further investigation should be followed to step forward to break the Kc values down to Ke and Kcb.
Cumulative stats can be summarized in a table.
Step 6. Calculate WUE using yd & the available data sets created
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 149.38722 <.0001
## sys 1 9 69.35126 <.0001
## spp 1 9 31.79966 0.0003
## sys:spp 1 9 0.07603 0.7890md.wue2 <- lme(data = wue, wue ~ sys + spp, random = ~1 | zn)
anova(md.wue2) # p < .001 for sys*** & .0001 for spp***## numDF denDF F-value p-value
## (Intercept) 1 10 149.38722 <.0001
## sys 1 10 76.41148 <.0001
## spp 1 10 35.03699 1e-04summary(md.wue2) ## Linear mixed-effects model fit by REML
## Data: wue
## AIC BIC logLik
## 72.61234 75.43708 -31.30617
##
## Random effects:
## Formula: ~1 | zn
## (Intercept) Residual
## StdDev: 1.228245 1.946488
##
## Fixed effects: wue ~ sys + spp
## Value Std.Error DF t-value p-value
## (Intercept) 16.710995 1.0428562 10 16.024257 0e+00
## sysOF -8.507480 0.9732438 10 -8.741366 0e+00
## sppTom -5.760829 0.9732438 10 -5.919205 1e-04
## Correlation:
## (Intr) sysOF
## sysOF -0.467
## sppTom -0.467 0.000
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.62774942 -0.46920278 0.06015411 0.27735282 1.60420132
##
## Number of Observations: 16
## Number of Groups: 4lsmeans(md.wue2, pairwise ~ sys, adjust = "tukey") # ht:of = a:b ## $lsmeans
## sys lsmean SE df lower.CL upper.CL
## HT 13.83 0.922 3 10.90 16.77
## OF 5.32 0.922 3 2.39 8.26
##
## 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 8.51 0.973 10 8.741 <.0001
##
## Results are averaged over the levels of: spplsmeans(md.wue2, pairwise ~ spp, adjust = "tukey") # pep:tom = a:b## $lsmeans
## spp lsmean SE df lower.CL upper.CL
## Pep 12.5 0.922 3 9.52 15.39
## Tom 6.7 0.922 3 3.76 9.63
##
## 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 5.76 0.973 10 5.919 0.0001
##
## Results are averaged over the levels of: syspwue <- 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 = .0098 for sys**## numDF denDF F-value p-value
## (Intercept) 1 3 145.78700 0.0012
## sys 1 3 34.69681 0.0098summary(md.p1) # regression model summary## Linear mixed-effects model fit by REML
## Data: pwue
## AIC BIC logLik
## 38.30871 37.47575 -15.15435
##
## Random effects:
## Formula: ~1 | zn
## (Intercept) Residual
## StdDev: 1.519099 1.974915
##
## Fixed effects: wue ~ sys
## Value Std.Error DF t-value p-value
## (Intercept) 16.570156 1.245788 3 13.30094 0.0009
## sysOF -8.225802 1.396476 3 -5.89040 0.0098
## Correlation:
## (Intr)
## sysOF -0.56
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.23291615 -0.47987147 0.04986552 0.40005615 1.25298429
##
## Number of Observations: 8
## Number of Groups: 4lsmeans(md.p1, pairwise ~ sys, adjust = "tukey") # ht:of = a:b## $lsmeans
## sys lsmean SE df lower.CL upper.CL
## HT 16.57 1.25 3 12.61 20.5
## OF 8.34 1.25 3 4.38 12.3
##
## d.f. method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## HT - OF 8.23 1.4 3 5.890 0.0098md.t1 <- lme(data = twue, wue ~ sys, random = ~1 | zn)
anova(md.t1) # ANOVA table, p = .0071 for sys** ## numDF denDF F-value p-value
## (Intercept) 1 3 56.23398 0.0049
## sys 1 3 43.50265 0.0071summary(md.t1) # regression model summary## Linear mixed-effects model fit by REML
## Data: twue
## AIC BIC logLik
## 37.16115 36.32818 -14.58057
##
## Random effects:
## Formula: ~1 | zn
## (Intercept) Residual
## StdDev: 1.189093 1.884536
##
## Fixed effects: wue ~ sys
## Value Std.Error DF t-value p-value
## (Intercept) 11.091005 1.114161 3 9.954580 0.0022
## sysOF -8.789159 1.332568 3 -6.595654 0.0071
## Correlation:
## (Intr)
## sysOF -0.598
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.98241024 -0.39980580 -0.03735849 0.29532909 1.61650006
##
## Number of Observations: 8
## Number of Groups: 4lsmeans(md.t1, pairwise ~ sys, adjust = "tukey") # ht:of = a:b## $lsmeans
## sys lsmean SE df lower.CL upper.CL
## HT 11.1 1.11 3 7.55 14.64
## OF 2.3 1.11 3 -1.24 5.85
##
## d.f. method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## HT - OF 8.79 1.33 3 6.596 0.0071To 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.
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.↩