Upscaling Evapotranspiration with Parsimonious Models in a North Carolina Vineyard

Abstract

Water stress can positively or negatively impact grape yield and yield quality, and there is a need for wine growers to accurately regulate water use. In a four-year study (2010–2013), energy balance fluxes were measured with an eddy-covariance (EC) system in a North Carolina vineyard (Vitis vinifera cv. Chardonnay), and evapotranspiration (ET) and the Crop Water Stress Index (CWSI) calculated. A multiple linear regression model was developed to upscale ET using air temperature (Ta), vapor pressure deficit (VPD), and Landsat-derived Land Surface Temperature (LST) and Enhanced Vegetation Index (EVI). Daily ET reached values of up to 7.7 mm day⁻¹, and the annual ET was 752 ± 59 mm, as measured with the EC system. The grapevine CWSI was between 0.53–0.85, which indicated moderate water stress levels. Median vineyard EVI was between 0.22 and 0.72, and the EVI range (max–min) within the vineyard was 0.18. The empirical models explained 75%–84% of the variation in ET, and all parameters had a positive linear relationship to ET. The Root Mean Square Error (RMSE) was 0.52–0.62 mm. This study presents easily applicable approaches to analyzing water dynamics and ET. This may help wine growers to cost-effectively quantify water use in vineyards.

1. Introduction

Quantifying water use in a vineyard is critical because plant available water in the root zone determines yield quantity and quality. The vines require no-stress conditions from bud break subsequent to fruit set, while moderate water stress is required during fruit ripening. Moderate water stress has an inhibitive effect on vegetative growth, which can enhance or positively impact berry composition, and eventually wine quality. Depending on the grapevine phenological stage, excessive water stress or oversupply of water can negatively impact yield and quality [1]. Water use in vineyards can be regulated by estimating the daily actual evapotranspiration (ET), i.e., the amount of water transpired by the crop, or evaporated from the soil and exposed surfaces. The goal in semi-arid wine regions is to reduce ET [2,3,4], however, wine growers in North Carolina (NC) and other humid regions seek to increase ET to a maximum because the high annual precipitation (P) inhibits vines to enter any water stress period during fruit ripening. A variety of techniques have been used to estimate ET in vineyards. Estimating latent heat flux (LE) using the Eddy covariance (EC) method has the advantage of generating continuous long-term ET datasets in vineyards on a sub-field or field level [5].

The Estimated Crop Water Stress Index (CWSI) is another method of quantifying plant water stress based on canopy temperature (Tc) measurements [6,7]. CWSI is related to the ratio of ET and potential ET (ET_p) [7]:

$$\mathrm{CWSI}=1-\frac{\mathrm{ET}}{{\mathrm{ET}}_{\mathrm{p}}}$$

(1)

The ET and Tc are related, because evaporative cooling of a transpiring plant reduces Tc, and water stressed plants close their stomates reducing ET, which leads to increased Tc. Several studies have successfully utilized the CWSI approach to quantify water stress [8,9]. More recently, the CWSI was applied in vineyards to map vine water stress and demonstrated that CWSI was related to stem water potential and yield parameters in grapevines [10,11,12,13].

There is a growing demand to cost-effectively quantify ET for whole vineyards. Previous studies used different modelling approaches or included remotely sensed data to estimate large-scale ET. One approach was the use of air-borne measured land surface temperature (LST) and air temperature (Ta) to estimate sensible heat (H). The LE and ET were then calculated as the residual energy in the energy balance equation [14]. Applying this concept, models for large-scale ET estimations have been developed. The Surface Energy Balance Algorithm for Land (SEBAL) estimates H using LST on very wet and dry conditions, i.e., at potential and zero LE [15]. A similar approach is the Mapping Evapotranspiration with Internalized Calibration (METRIC), which combines the SEBAL approach with micro-meteorological measurements on the ground [16]. The METRIC model was used to estimate regional vineyard ET in Chile in comparison with EC measurements [2]. The Simplified Surface Energy Balance Index (S-SEBI) model [17], another energy balance model based on LST, had previously been utilized to simulate ET in a vineyard in the Mediterranean region [4]. The Alexi/Disalexi approach uses a two-source energy balance (TSEB) model together with satellite-derived Normalized Difference Vegetation Index (NDVI), LST, and other meteorological data in a data fusion method to create daily ET products [3,18]. This approach was used for tempo-spatial estimation of ET in different agroecosystems, such as in vineyards in California [3,18]. More methods of ET estimation have been described in the literature [19].

While these modeling approaches are well developed to estimate ET, a more parsimonious model to upscale ET from sub-field ground measurements to total production areas is more accessible to growers. Volumetric soil water content (VWC) as a predictor of ET, and measured with EC stations, was used to upscale ET in the North American monsoon region [20]. Another way to estimate ET was the use of the reference evapotranspiration (ET₀) multiplied with a crop coefficient (kc) [21], of which the former was calculated from meteorological data, and the latter was substituted by a satellite-derived vegetation index (VI) [14,22]. In addition, [23] reported that substituting ET₀ with Ta showed good agreement in upscaling ET, using a simple empirical model with the Enhanced Vegetation Index (EVI) as a kc substitute in riparian vegetation in the Southwest U.S. Other climatic variables may also be used to substitute ET₀. Landsat imagery provides EVI and LST products at a 30 × 30 m resolution, which would allow ET to be estimated on a subfield level using parsimonious regression models.

In this four-year study, actual ET and CWSI were estimated in a commercial vineyard in NC, USA, using the EC method as a comparison to quantify associated ET. This study’s aim was to monitor the water dynamics and quantify ET, and to eventually upscale ET from sub-field to vineyard level by developing a parsimonious empirical model using micro-meteorological data and Landsat 7 imagery.

2. Materials and Methods

2.1. Study Site

The study was conducted from April 2010 to November 2013 at a commercial vineyard (~8.7 ha) near Dobson, NC (−80.777° W; 36.357° N; 366 m a.s.l.). The soil type was a fine, kaolinitic, mesic, Typic Kanhapludult, and the soil texture was sandy clay loam. The vines (Vitis vinifera cv. Chardonnay) were planted in 2001 in the north to south direction at 1.8 m spacing in 2.7 m wide rows. The vines were cordon-trained and spur-pruned. The canopy height was 0.9 to 1.9 m from the ground. Canopy width was 0.3 to 0.8 m. The vine rows were periodically treated with herbicide to control weeds under vine vegetation. The interrows were seeded with tall fescue (Festuca arundinacea Shreb), yet native vegetation was prevalent [24]. The vineyard was not irrigated. The climate was humid with mean daily Ta ranging from −2.2 to 27.0 °C (Figure 1A) and mean annual precipitation (P) of 927 mm (Figure 1G).

2.2. Eddy Covariance System

One EC station was installed on the southern side of the vineyard from April 2010 to November 2013, and turbulent fluxes of LE and sensible heat flux (H) were computed using the EC method. Fetch-to-height ratio was ≥40:1 in all directions. The EC station was equipped with a high frequency open-path infrared gas analyzer (IRGA), that measured water vapor (LI-7500, LICOR Biosciences Inc., Lincoln, NE, USA). High frequency instantaneous wind speed velocity components and sonic temperature were measured with a 3-D sonic anemometer (CSAT, Campbell Scientific, Inc., Logan, UT, USA). The instantaneous data were collected at 20 Hz, and 15-min averages were generated. The H and LE fluxes were computed as the covariance between the instantaneous vertical wind velocity, sonic temperature, and water vapor density. The LE and H fluxes were coordinate rotated [25], and LE was additionally corrected for variations in temperature (Ta) and air density based on the Webb–Pearman–Leuning algorithm [26]. Ancillary instrumentation included a net radiometer to measure net radiation (Rn; Q7 net radiometer, Radiation Energy Balance, Seattle, WA, USA), an air temperature and relative humidity (rH) probe (Vaisala HMP45c, Vaisala, Vantaa, Finland), a tipping bucket rain gauge (Model TR525 USW, Texas Electronics, Dallas, TX, USA) to measure P, soil moisture probes (CS615, Campbell Sci., Logan, UT, USA) at 6 cm soil depth to measure VWC, thermocouples to measure soil temperature, heat flux plates (HFT-1, Radiation Energy Balance, Seattle, WA, USA) to measure soil heat flux (G), and four infrared temperature (IRT) sensors (IRTS-P, Apogee Instruments Inc., Logan, UT, USA) above the vines and fescue at 1.5 m, and above the grapevine row and grass interrows in approximately 6.1 m to measure Tc (referred to as Tc_v, Tc_f, Tc_r, Tc_ir, respectively).

The datasets were screened for sensor failures and statistical outliers and subsequently interpolated to fill data gaps. All data that violated empirically set upper and lower thresholds were excluded from further analysis. Additional screening features of the turbulent fluxes included periods of precipitation [27], possible water condensation on the instruments [28], low wind turbulence, wind direction opposite to the sensor head direction [29], and the manufacturer’s internal IRGA and CSAT warning system. We also used the water vapor saturation measurement from the humidity probe to screen LE, since the humidity sensor is not affected by dust or rainfall events [30]. Data were also screened using the interquartile range (IQR) procedure, where outliers were defined as 1st/3rd Quartile ± 3.5 × IQR. The 15-minute dataset was interpolated (or gap filled) following an inverse weighted time average procedure [29]. Daily values were calculated as the daily average multiplied by the number of samples per day. The Bowen ratio (β) was calculated as the quotient of H and LE. The ET was calculated as the quotient of LE and the latent heat of vaporization (LV), of which the latter was calculated as proposed by [31]. Further information on outlier screening and interpolation procedures can be found in [32]. The annual ET was calculated as mean daily ET multiplied by the number of days year⁻¹.

2.3. Crop Water Stress Index (CWSI)

The CWSI was calculated as [6,7]:

$$\mathrm{CWSI}=\frac{\mathrm{dT}-{\mathrm{dT}}_{\min }}{{\mathrm{dT}}_{\max }-{\mathrm{dT}}_{\min }}$$

(2)

Where: CWSI = Crop Water Stress Index of vines, fescue, grapevine row, and interrow (referred to as CWSI_v, CWSI_f, CWSI_r, CWSI_ir, respectively); dT = the difference between Tc_v, Tc_f, Tc_r, Tc_ir and Ta (°C); dT_min = the lower, water-stress free limit (°C), dT_max = the upper, highly water-stressed limit (°C).

Data from full sun middays (11 AM–3 PM; R_n > 500 W m⁻²) between May and October were used to calculate CWSI, where vines had full foliage, and cloud cover did not influence Tc. The boundaries dT_min and dT_max describe the crop- and location-specific limits of non-stressed (plant at potential transpiration) and fully water stressed (virtually no transpiration) conditions of a crop canopy. The lower boundary dT_min was calculated as [7,33]:

$$\begin{array}{ll}d{T}_{min}& =\frac{r_{ae}\times \left(R_n-G-S_G\right)}{\left(\rho \times c_p\right)}\\ & \times \frac{\gamma \times \left(1+\frac{r_{cp}}{r_{ae}}\right)}{\Delta +\gamma \times \left(1+\frac{r_{cp}}{r_{ae}}\right)}-\frac{\left(e_s-e_a\right)}{\Delta +\gamma \times \left(1+\frac{r_{cp}}{r_{ae}}\right)}\end{array}$$

(3)

where: r_ae = effective aerodynamic resistance (s m⁻¹), R_n = net radiation (MJ m⁻² s⁻¹), G = soil heat flux (MJ m⁻² s⁻¹), S_G = soil heat storage (MJ m⁻² s⁻¹), ρ = air density (kg m⁻³), c_p = heat capacity of air (MJ kg⁻¹ °C⁻¹), r_cp = crop canopy resistance at potential transpiration (s m⁻¹), Δ = slope of the saturated vapor pressure-temperature relation (kPa °C⁻¹), γ = psychrometric constant (kPa °C⁻¹), e_s = saturated vapor pressure (kPa); e_a = actual vapor pressure (kPa).

The upper boundary dT_max was calculated as [7,33]:

$$d{T}_{max}=\frac{r_{ae}\times \left(R_n-G-S_G\right)}{\left(\rho \times c_p\right)}$$

(4)

e_s was calculated as a function of Ta [34], c_p and ρ were calculated as described in [35], and γ as described in [21]. The e_a was calculated as e_s multiplied with rH. The Δ was calculated with the average of Tc (that is, Tc_v, Tc_f, Tc_r, Tc_ir) and Ta [33]. The r_cp was set to 5 s m⁻¹ [7], and r_ae was calculated according to [36]. The S_G was calculated using the calorimetric method [37]. The CWSI values were screened for outliers outside its limits of 0 (no water stress) and 1 (maximum water stress), and subsequently gapfilled using the inverse weighted time average procedure.

2.4. Enhanced Vegetation Index (EVI) and Land Surface Temperature (LST) from Landsat 7

The EVI and LST Level-2 product of Landsat 7 satellite images were provided by the U.S. Geological Survey (USGS) [38,39]. The Landsat 7 EVI product had a resolution of 30 × 30 m. Landsat 7 LST was calculated with the thermal band which had a 60 × 60 m resolution. The final LST Level-2 product was resampled to a 30 × 30 m resolution. Note that USGS Landsat Surface Temperature Science Product may have reported unvalidated results for certain observational conditions [40]. In total, 37 EVI and 32 LST scenes with little or no cloud cover between May 2010 to October 2013 were requested. Five LST products were not available, probably due to the current provisional nature of LST. The EVI product as provided by USGS was calculated as:

EVI = 2.5 × ((ρ_NIR − ρ_Red)/(ρ_NIR + 6 × ρ_Red − 7.5 × ρ_Blue + 1))

(5)

Where: EVI = Enhanced Vegetation Index, ρ_NIR, ρ_Red, and ρ_Blue were the surface reflectance of the near-infrared, red, and blue wavelength bands, respectively.

The Scan Line Corrector (SLC) of the Landsat 7 failed in 2003, which caused an unequal ground track of the satellite with duplicate and missing imaged areas. Landsat scenes with missing data were not included in this study.

The EVI and LST scenes were cropped to field boundaries of the vineyard as well as to the 80% daytime footprint area (i.e., two to three data points), assuming that the greatest share of LE was derived from this footprint. The median EVI and LST around the EC station was then calculated and used for further analysis.

2.5. Statistical Analysis and Data Preparation

In this study a multiple linear regression analysis with on-site Ta and VPD, as well as LST and EVI within the EC footprint as the independent and ET as the dependent variable was established (p < 0.05). These empirical and parsimonious models were then used to upscale ET for the whole vineyard with vineyard EVI and on-site Ta, VPD or vineyard LST, respectively (referred to as ET_Ta, ET_LST, and ET_VPD, respectively). The relationship between ET and CWSI was analyzed with linear regression analysis. The assumption was that Ta and VPD measured at the EC station did not differ substantially within the vineyard. All data preparation, graphing, and analysis was carried out using R [41].

3. Results

Daily ET, Ta, and VPD showed the typical seasonal pattern with high values during summer and low values during winter (Figure 1A,B,D). The VPD reached values of up to 1.90 kPa. Maximum ET in the growing season (i.e., May to October) ranged from 5.6 mm day⁻¹ in 2013 to 7.7 mm day⁻¹ in 2012. During the growing season median β was between 0.36 and 0.40, and in the off-season, it was between 0.47–0.84 (Figure 1C). The VWC in the surface 6 cm of soil under tall fescue ranged from 0.29 m³ m⁻³ in 2010 and 0.42 m³ m⁻³ in 2013, and under grapevines VWC was between 0.24–0.35 m³ m⁻³. The VWC pattern closely followed rainfall-events (Figure 1E–G). Mean annual ET was 752 ± 59 mm year⁻¹. Note that the small annual ET in 2010 was related to the late onset of the study and a data gap from DOY 153–175 (

Table 1. Annual evapotranspiration (ET) from 2010–2013.

Year	ET (mm year−1)
2010	579
2011	837
2012	807
2013	784
Mean	752 ± 59

, Figure 1B). Precipitation (P) was above ET in 2010 –2013 with a mean P-ET ratio of 1.21.

The dTmin values of fescue, grapevines, and measurements at 6.1 m above row and interrow was significantly related to VPD (p < 0.05), explaining 93–95% of the variation. The median CWSI_f was lowest in 2011 with 0.81 and highest in 2010 with 1.0, and for CWSI_v was lowest in 2010 with 0.53 and highest in 2013 with 0.85, respectively. The annual median CWSI_ir and CWSI_r ranged between 0.76 and 0.91 and between 0.77 and 0.91, respectively. The CWSI was highly variable among years (Figure 2). There was a significantly negative linear relationship between ET and CWSI_r (p < 0.05), i.e., with increasing CWSI the ET decreased. However, the relationship could only explain 15% of the variation.

Vineyard median EVI ranged from 0.22–0.72 (Figure 3), with a mean within-vineyard range (calculated as the difference of maximum and minimum EVI) of 0.18. Median vineyard LST ranged from 16.0–40.1 °C from 2010 to 2013, with a mean range within the vineyard of 4.44 °C (Figure 4). Median EVI within the EC footprint followed a seasonal pattern and ranged from 0.28–0.65 (Figure 5A). Median LST within the EC footprint also followed a seasonal pattern and ranged from 0.8 to 42.6 °C (Figure 5B). The LST was significantly related to Tc at overpass time and lowest RMSE of 1.13 °C was measured with an IRT sensor at 6.1 m above the grapevine row (

Table 2. The RMSE (°C), and R² of linear regression models with Tc (°C) measured with IRT sensors at time of satellite overpass as independent variable, and LST (°C) as dependent variable.

Tc	R2	RMSE (°C)
Tcf	0.99	2.18
Tcv	1.00	1.68
Tcr	1.00	1.13
Tcir	1.00	1.49

IRT = infrared temperature; RMSE = root mean squared error, R² = coefficient of determination; LST = Land surface temperature, Tc = canopy temperature.

).

Footprint EVI and on-site Ta, LST, and VPD were used as independent variables to establish multiple linear regression models with ET as the dependent variable (

Table 3. Coefficients of a multiple linear regression analysis with EVI and Ta, VPD, or LST as independent and ET as dependent variable. The RMSE refers to ET measured with the EC station versus modeled ET within the 80% footprint area.

Model	Coefficients	Mean	SE	p	Adj. R2	dF	Variable range	Adj. R2 Measured vs. Modeled ET	RMSE Measured vs. Modeled ET
ETTa	Intercept	−1.629	0.592	<0.05	0.76	28		0.77	0.60
	EVI	7.181	1.821	<0.001			0.28–0.65
	Ta (°C)	0.069	0.025	<0.01			−4.55–29.32
ETVPD	Intercept	−2.131	0.457	<0.001	0.84	28		0.84	0.52
	EVI	7.866	1.159	<0.001			0.28–0.65
	VPD (kPa)	1.645	0.332	<0.001			0.23–1.89
ETLST	Intercept	−2.438	0.631	<0.001	0.75	24		0.76	0.62
	EVI	8.877	2.017	<0.001			0.28–0.63
	LST (°C)	0.042	0.021	<0.1			2.81–40.15

Adj. R² = adjusted coefficient of determination, dF = degrees of freedom, EVI = Enhanced Vegetation Index, LST = Land surface temperature, RMSE = root mean squared error, Ta = air temperature, SE= standard error, p = significance level, VPD = vapor pressure deficit.

). All models were significant at p < 0.05, explaining 75% to 84% of the variation. All parameters had a positive linear relationship to ET, i.e., with increasing EVI, Ta, LST or VPD, the daily ET increased and vice versa. A linear regression analysis between calculated ET and measured ET was significant at p < 0.05 and had a root mean square error (RMSE) of 0.52–0.62 mm (Table 3). Note that combining more independent variables did not improve the empirical models, so the more parsimonious models were chosen. While the overall EVI-LST model was significant, the LST parameter was only significant at p < 0.1, due to higher variation in the LST dataset. The derived coefficients were then used to calculate vineyard ET. Median daily vineyard ET_Ta, ET_LST, and ET_VPD ranged from 0.4–4.8, 0.7–4.6, and 0.8–5.1 mm day⁻¹, respectively (Figure 5C–E). The mean within-vineyard ET range was 1.31, 1.67, and 1.48 mm for ET_Ta, ET_LST, and ET_VPD, respectively (Figure 6).

4. Discussion

In this study, ET was measured using an EC system in a NC vineyard. Similar daily ET values for interrrow areas had previously been measured with a micro-Bowen Ratio energy balance system (MBR) [24], showing that the estimated ET values were reasonable. The authors also showed that the interrow tall fescue increased ET compared to bare soil, and that the contribution of tall fescue to ET was higher in the off season than during the growing season [24]. Note that ET measured with the EC system was a shared signal of both, vines and tall fescue, as the EC footprint was greater than for the MBRs.

The Landsat 7 EVI and LST were as well based on the reflectance of both, vines and grass at a resolution of 30 × 30 m. The LST was best correlated with Tc_r and Tc_ir (Table 2), probably owing to the higher IRT footprint measuring a combined signal of interrow, soil, and vines. Significant relationships with Landsat brightness temperature and IRT sensor measurements [42], and Landsat LST and soil temperature [43] had previously been reported. The ET values presented were specific to the local environment as well as vineyard management (as described by LST, EVI, and EC derived data), and any upscaling approaches for other systems may need additional adjustments.

The VWC at 6 cm soil depth was mostly above the wilting point of 0.27 m³ m⁻³ during the study period [44], especially under tall fescue, which would indicate no severe water-stress soil conditions. Yet, CWSI values were relatively high, an indicator that water-stress occurred, at least during high-radiation time periods during midday. In a Mediterranean vineyard CWSI values from 0.45−0.65 were reported, which was lower than in this study [12]. The authors also found a positive relation between yield parameters and CWSI, the higher the CWSI the higher the berry weight. Assuming a similar yield response, mid-season CWSI would indicate water stress levels that would benefit yield but would be too high at the onset of the growing season, when water stress is not beneficial. Yet, yield also depends on berry and cluster number, and other factors such as trimming can further influence yield. Thus, further analysis with a larger dataset including yield data would be needed to confirm this relationship for our study. The linear relationship between ET and CWSI_v would indicate the effect of evaporative cooling on Tc (Equation 1), and similar relationships have previously been reported for other agro-ecosystems [45]. Yet, CWSI was a weak predictor of ET compared to EVI, Ta, LST, or VPD in this study. The footprint of EC derived ET was greater than IRT-measured Tc. In addition, the CWSI was calculated for midday time periods on days with high net radiation, that is periods of maximum water stress. This may explain the weak relationship with EC-measured daily ET values.

In a previous study, EVI and Ta were used as substitutes for kc and ET₀ to estimate ET from riparian vegetation [14]. In this study, ET_Ta showed good agreement with measured ET, demonstrating that EVI and Ta can be used to estimate ET. However, the best ET estimates were achieved with VPD as the independent variable. [14] found a sigmoidal relationship with Ta and an exponential relationship with EVI. This was not the case in this study, where a linear relationship was used. There was no upper asymptote in the ET-Ta relationship, i.e., there was still an ET response even at maximum Ta. However, this may not be the case for exceptionally dry and hot years, where the model presented may overestimate ET.

Estimating ET_Ta requires only a point measurement of Ta and EVI maps. Yet, higher accuracy (i.e., lowest RMSE) was achieved with VPD as the independent variable, however, this would require two measurements, Ta and rH or water vapor density. ET_Ta and ET_VPD estimation assumed that Ta and rH did not substantially change within the vineyard. However, even on smaller scales (field-level) we found substantial variation in LST, and it may be recommended to use satellite-derived LST for ET upscaling [3,4]. The drawback using Landsat VI and LST products was that Landsat had a 16-day overpass-period, and cloud cover could impact the image quality. Data fusion approaches or other interpolation methods may then be needed to generate datasets in shorter time periods [3]. Despite these caveats, it was possible to accurately upscale ET for the whole vineyard with these parsimonious models.