Evapotranspiration Estimate over an Almond Orchard Using Landsat Satellite Observations

Abstract

California growers face challenges with water shortages and there is a strong need to use the least amount of water while optimizing yield. Timely information on evapotranspiration (ET), a dominant component of crop consumptive water use, is critical for growers to tailor irrigation management based on in-field spatial variability and in-season variations. We evaluated the performance of a remote sensing-based approach, Mapping Evapotranspiration at high Resolution with Internalized Calibration (METRIC), in mapping ET over an almond orchard in California, driven by Landsat satellite observations. Reference ET from a network of weather stations over well-watered grass (ET_o) was used for the internal calibration and for deriving ET at daily and extended time period, instead of alfalfa based reference evapotranspiration (ET_r). Our study showed that METRIC daily ET estimates during Landsat overpass dates agreed well with the field measurements. During 2009–2012, a root mean square error (RMSE) of 0.53 mm/day and a coefficient of determination (R²) of 0.87 were found between METRIC versus observed daily ET. Monthly ET estimates had a higher accuracy, with a RMSE of 12.08 mm/month, a R² of 0.90, and a relatively small relative mean difference (RMD) of 9.68% during 2009–2012 growing seasons. Net radiation and Normalized Difference Vegetation Index (NDVI) from remote sensing observations were highly correlated with spatial and temporal ET estimates. An empirical model was developed to estimate daily ET using NDVI, net radiation (R_n), and vapor pressure deficit (VPD). The validation showed that the accuracy of this easy-to-use empirical method was slightly lower than that of METRIC but still reasonable, with a RMSE of 0.71 mm/day when compared to ground measurements. The remote sensing based ET estimate will support a variety of State and local interests in water use and irrigation management, for both planning and regulatory/compliance purposes, and it provides the farmers observation-based guidance for site-specific and time-sensitive irrigation management.

1. Introduction

Almond has become the first leading revenue-generating crop commodity in California, with a total value of annual production over 5.9 billion dollars in 2014. California’s almond acreage has been rapidly growing, from 700,000 acres (bearing acreage: 590,000 acres) in 2005 to 1,020,000 acres (bearing acreage: 870,000 acres) in 2014, according to the California Department of Food and Agriculture (CDFA) [1]. On the other hand, almond growers face water supply challenges due to the increasing frequency and severity of recurring droughts. Water use by almond orchards has not been well studied, especially under various cropping conditions, and this knowledge gap limits the capacity for growers to achieve resource-efficient water management.

Agricultural consumptive water use consists of water removed by evapotranspiration (ET), via soil evaporation (E) and plant transpiration (T). Accurate estimation and mapping of crop ET across space and time are essential for improving field-scale irrigation management and for watershed and regional water planning and management. Ground-based flux measurement methods such as eddy covariance (EC) [2], Bowen ratio (BR) [3], and most recently surface renewal (SR) [4], can provide half-hourly, hourly, and daily ET measurements and thus offer guidance to growers for time sensitive irrigation scheduling. However, these measurements are limited by the small scale of the footprint area of the measurement stations and they cannot capture heterogeneity within each agriculture field, and sometimes cannot represent the entire field or orchard conditions. The equipment for measuring field ET through ground methods is usually expensive, while their installation and maintenance could also be labor-intensive.

In agricultural production areas, ET has been traditionally estimated as a product of the weather-based ET of a hypothetical reference surface, such as alfalfa (ET_r) or grass (ET_o) under well-watered condition, and crop coefficient (K_c) [5,6]. K_c is often generalized for each crop type and over a few growth stages from limited field measurements, and thus does not capture the detailed temporal dynamics of crop growth and spatial variability of crops within and across the field blocks. For example, values of almond K_c of 0.40 at initial stage, 0.90 at mid-season and 0.65 at the end of the season are usually considered for irrigation management in all almond orchards [5]. Near real time values of ET_o are readily accessible to growers throughout the State of California from the California Irrigation Management Information System (CIMIS) network of more than 150 weather stations, and they are used extensively for irrigation management in California. However, spatially-explicit seasonal K_c values for estimating water use by almond orchards are not yet available in the literature. Updating the almond water use information is particularly important for California’s almond growers to improve their on-farm irrigation management, since almond production industry experienced major changes in cultivars, rootstocks, planting density, canopy management practices, and irrigation technology and methods during the past 15 to 20 years.

New approaches to take advantage of remote sensing observations to map ET have been developed, especially in recent decades, due to increasingly available free satellite observations from multiple sensors [7,8,9,10]. Empirical algorithms were built to relate ET to remotely sensed vegetation indices and weather data [11,12,13,14]. These algorithms typically perform well in areas where they are calibrated, but the uncertainty usually increases when applied to other areas. The semi-empirical Priestley–Taylor equation, calibrated with ET measurements, has also shown promising results for estimating ET at continental to global scale [7,8]. A suite of algorithms that were based on surface residual of the energy balance (REB) data, using both optical and thermal remote sensing observations, has also provided estimates of ET. One-source REB models include Surface Energy Balance Algorithm for Land (SEBAL) [15], Surface Energy Balance System (SEBS) [16], Mapping Evapotranspiration at high Resolution with Internalized Calibration (METRIC) [17], and Simplified Surface Energy Balance (SSEB) [18,19]. Two-source models, such as the Two-Source Model (TSM) [20,21] and the Disaggregated Atmosphere-Land Exchange Inverse (DisALEXI) [22,23,24,25], integrate the biophysical processes and remote sensing data to estimate the plant transpiration and evaporation from non-vegetated surfaces separately.

METRIC estimates the near surface temperature gradient (dT) as an indexed function of radiometric surface temperature, and thus does not rely on accurate estimate of land surface temperature from thermal imagery [17,26,27]. It also uses hourly reference ET to auto-calibrate the sensible heat calculation for each satellite image. This internal calibration makes ET estimates more robust. METRIC has been widely used to estimate agricultural water use in Idaho, California, and New Mexico for management of water rights, irrigation scheduling, and water resource planning [28,29]. Regional scale ET mapping has also been done in the Texas High plains [30], in the mountainous areas of northern Portugal [31], and in the oasis area in Heibe River Basin in China [32], and the results showed reasonable accuracy. The METRIC method was used to provide good estimates of ET over olive orchards [33,34,35] and over cotton fields [36]. It tends to outperform TSEB when sufficient ancillary data are not available [36]. To our knowledge, however, there are no studies on almond ET estimation by METRIC, and it is not clear how remote sensing-based ET approaches perform in almond orchard with the complex canopy structure of almond trees.

In this paper, we aim to: (1) assess the accuracy of ET estimates at the field scale in an almond orchard using the METRIC approach, driven by Landsat satellite observations and ET_o from CIMIS; (2) examine the dominant drivers for temporal and spatial variation in crop water use; and (3) develop and assess empirical methods that are easy to use for ET estimation at the farm scale.

2. Study Area and Datasets

2.1. Study Area

The study area was a 60 ha mature almond orchard (Prunusdulcis) located near Lost Hill in Kern County, in the southern San Joaquin Valley of California (35.510°N, 119.667°W; Figure 1). The almond trees were planted in 1999 on a Milham sandy loam, with South to North orientation, 6.4 m between trees in the rows, and 7.3 m between the rows [37]. Two types of micro-irrigation systems, i.e., drip and fanjet, were used to water the orchard. This orchard was relatively homogeneous in terms of canopy structure since similar soil water storage was achieved for all trees on purpose through a careful management of the two micro-irrigation systems during 2009–2012. Almond trees were about 7 m tall during our study period from 2009 to 2012.

San Joaquin Valley has a typical Mediterranean climate. Mean annual temperatures ranged from 15.9 to 17.2 °C during 2009–2012. Summer months are hot and dry, while winter months are mild and wet (Figure 2a). Most rainfall occurs in winter and spring from November to April. Precipitation varied significant from year to year, e.g., annual precipitation varied from 98.5 mm in 2009 to 235.7 mm in 2010. Almond trees flower in late February to early March and reach almost full canopy cover by late March. Harvest typically occurs in late August to mid-September.

The reference ET_o related to grass had a strong seasonal cycle, ranging from 23.8 mm/month in January to 214.1 mm/month in July (Figure 2b). The annual accumulated reference ET_o was 1492.2, 1393.8, 1330.7, and 1450.8 mm/year in 2009–2012, respectively. The almond orchard was irrigated following typical practices in the upper San Joaquin Valley.

2.2. Remotely Sensed Data

Landsat 5 TM and 7 ETM+ observations were obtained from USGS (http://earthexplorer.usgs.gov/) for path 42 and row 35 during 2009–2012, covering the almond orchard study area (Supplementary Materials Table S1). Forty-six cloud-free images were obtained during the four year study period, with an average of 12 images per year, mostly from April to November. The Landsat revisiting time was 16 days, but most images during December to March rain season were not usable for our analysis due to frequent cloud cover. All the images were subsetted to a common domain, as shown in Figure 1a. A careful examination of individual Landsat 7 ETM+ images showed that this domain had no data gaps during the study period. This subset of 801 by 801 pixels included both bare soil and well-watered full cover vegetation areas, as needed by METRIC approach for the purpose of selecting cold and hot pixels to calculate sensible heat.

2.3. Micrometeorological Measurements

2.3.1. Station Reference Evapotranspiration

Reference ET_o, associated with well-irrigated grass, and other meteorological data were obtained from CIMIS (http://www.cimis.water.ca.gov/), as ancillary inputs to the METRIC algorithm. CIMIS was jointly developed by the California Department of Water Resources (DWR) and the University of California, Davis (UCD) in 1982 to inform water resource and irrigation management. CIMIS currently operates a network of 164 automated stations measuring weather parameters over well-watered standardized grass surface across California. A typical ET_o station consists of several sensors including pyranometer, thermistor, HMP35 temperature and humidity probe, wind vane, anemometer, and tipping-bucket rain gauge to collect information on solar radiation, soil temperature, air temperature, relative humidity, wind direction and speed, and precipitation, respectively. Data are collected every minute and archived as totals or averages each hour. Hourly and daily values of reference ET_o are not measured directly, but calculated from meteorological data using both CIMIS Penman equation and the standardized reference evapotranspiration equation for short canopies [6]. The standardized ET_o equation data were used in this study. ET_o was only used for the internal calibration of energy balance and for calculating the fraction of grass reference ET as an index to represent the relative change in weather [17] (See 3.2 below).

The closest CIMIS station, Belridge station (CIMIS #146; Lat/Lon: 35.50°N/119.69°W), is approximately 2 km away from the almond orchard (Figure 1a). Belridge station was located near the center of footprint of Landsat 5 TM and 7 ETM+ image subsets, and thus was selected to provide weather input data sets and reference ET_o for the METRIC approach in the study area. In particular, we used hourly reference ET_o, air temperature, vapor pressure and wind speed at Landsat overpass time, typically around 10:30 a.m., and daily reference ET_o during Landsat overpass dates from Belridge CIMIS station. Hourly reference ET_o was used to internally calibrate METRIC approach for sensible heat flux calculation. Air temperature was required to correct land surface temperature (LST) estimates from Landsat thermal imagery. Daily reference ET_o was used to convert instantaneous ET at Landsat overpass time to daily ET over the almond orchard in the METRIC approach, as shown below.

2.3.2. Actual Evapotranspiration Measurements

Evapotranspiration was measured with a micrometeorological tower installed 9 m above the ground in the almond orchard in early 2008. This system includes a net radiometer, sonic anemometer, two thermocouples, and soil heat flux plates, which were assembled to measure components of surface energy fluxes [38]. A REBS Inc. Q7.2 net radiometer was mounted at about 10 m height on the flux tower to measure net radiation (R_n). Sensible heat flux (H) was measured with both an RM Young 81000RE 3-d sonic anemometer and with Campbell Scientific Inc. FW3 fine-wire thermocouples as described in Shapland et al. [39]. The ground heat flux data were measured using three sets of ground heat flux packages, including a heat flux plate (REBS Inc., Seattle, WA, USA, HFT3) and a soil temperature averaging sensor (Campbell Scientific Inc., Logan, UT, USA, Tcav), which were established near the tower to measure the average ground heat flux (G) at the soil surface [40]. Latent heat flux (LE) was then estimated as the residual of the energy balance (LE = R_n − G − H) every half hour and converted to ET by dividing the LE in MJ m⁻² per half hour by 2.45 MJ mm⁻¹ of water vaporized. Hourly and daily ET of the almond orchard was calculated from the half-hourly data collected from 2009 through 2012. Using the station and measurement methodology described in Shapland et al. [39], there were almost no missing half hour estimates of almond ET. Daily and monthly calculations of ET were used to validate the METRIC approach for estimating and mapping almond ET.

3. Methodology

3.1. Landsat 5 TM and 7 ETM+ Data Preprocessing

The digital number (DN) of visible, near infrared and shortwave infrared bands of Landsat 5 TM and 7 ETM+ was first converted to radiance, and then atmospherically corrected to surface reflectance by using Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [41]. Surface albedo [42], normalized difference vegetation index (NDVI) [43] and leaf area index (LAI) [44] were derived from the spectral reflectances to calculate energy fluxes (R_n, G, and H) in the METRIC approach.

The original resolution of the thermal band (band 6) is 120 m. We here downscaled the thermal data to a 30-meter resolution to match the resolution of the reflective bands by using cubic convolution method. The thermal band DN was then converted to radiance using the coefficients provided in the imagery metadata. A correction to the radiance was further performed using the method of Wukelic et al. [45], and the Plank equation was then applied to compute the land surface temperature (LST). Cloud fraction (CF) mask for each Landsat image, provided by USGS (http://earthexplorer.usgs.gov/), was also used to remove data affected by cloud and cloud shadow.

3.2. METRIC Evapotranspiration Estimate: General Approach

METRIC™, a satellite-based surface energy balance approach, has been widely used to map field-scale crop ET, mostly over row crops. It is a one-source model that considers transpiration of vegetation and evaporation of soil as a whole. Latent heat flux (LE), directly related to ET, is computed as a residual according to surface energy balance equation (Equation (1)):

$$\mathrm{LE}={\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-\mathrm{H}$$

(1)

The major energy components are estimated as below:

$${\mathrm{R}}_{\mathrm{n}}=\left(1-\mathsf{\alpha }\right)*{\mathrm{R}}_{\mathrm{s}}\downarrow +{\mathrm{R}}_{\mathrm{l}}\downarrow -{\mathrm{R}}_{\mathrm{l}}\uparrow -\left(1-{\mathsf{\epsilon }}_0\right){*\mathrm{R}}_{\mathrm{l}}\downarrow$$

(2)

$$\mathrm{G}=\left(0.05+0.18{*\mathrm{e}}^{-0.52*\mathrm{LAI}}\right){*\mathrm{R}}_{\mathrm{n}}\left(\mathrm{LAI}\ge 0.5\right)$$

$$\mathrm{G}=\left(1.80*\frac{\mathrm{Ts}-273.16}{{\mathrm{R}}_{\mathrm{n}}}+0.084\right)*{\mathrm{R}}_{\mathrm{n}}(\mathrm{LAI}<0.5)$$

(3)

$$\mathrm{H}=\left({\mathsf{\rho }*\mathrm{c}}_{\mathrm{p}}*\mathrm{dT}\right)/{\mathrm{r}}_{\mathrm{ah}}$$

(4)

where $\mathsf{\alpha }$ is the surface albedo, R_s↓ is the incoming shortwave radiation (W/m²), R_l↓ is the incoming longwave radiation (W/m²), R_l↑ is the outgoing longwave radiation (W/m²), ${\mathsf{\epsilon }}_0$ is the broad band surface emissivity, $\mathsf{\rho }$ is the air density (1.15 kg/m³), ${\mathrm{c}}_{\mathrm{p}}$ is the air specific heat (1004 J/kg/K). Ts is the land surface temperature (K), and dT is the near surface-air temperature difference (K). dT is assumed to have a linear relationship with Ts (Equation (5)) as developed by Bastiaanssen [46] (Figure S1).

$$\mathrm{dT}=\mathrm{a}*\mathrm{Ts}+\mathrm{b}$$

(5)

where the coefficients a and b were calculated using two sets of “hot” and “cold” anchor pixels, as shown below:

$$\mathrm{a}=\frac{{\mathrm{dT}}_{\mathrm{hot}}-{\mathrm{dT}}_{\mathrm{cold}}}{{\mathrm{Ts}}_{\mathrm{hot}}-{\mathrm{Ts}}_{\mathrm{cold}}}\mathrm{b}={\mathrm{dT}}_{\mathrm{hot}}-\mathrm{a}*{\mathrm{Ts}}_{\mathrm{hot}}$$

(6)

where ${\mathrm{dT}}_{\mathrm{hot}}=\frac{{\mathrm{H}}_{\mathrm{hot}}*{{\mathrm{r}}_{\mathrm{ah}}}_{\mathrm{hot}}}{\mathsf{\rho }*{\mathrm{c}}_{\mathrm{p}}}$ and ${\mathrm{dT}}_{\mathrm{cold}}=\frac{{\mathrm{H}}_{\mathrm{cold}}*{{\mathrm{r}}_{\mathrm{ah}}}_{\mathrm{cold}}}{\mathsf{\rho }*{\mathrm{c}}_{\mathrm{p}}}$. The initial values of H, ${\mathrm{r}}_{\mathrm{ah}}$, and Ts were from the selected cold and hot pixels at the beginning of the iteration.

${\mathrm{r}}_{\mathrm{ah}}$ is the aerodynamic resistance to heat transport (s/m). An iterative method was applied to calculate both ${\mathrm{r}}_{\mathrm{ah}}$ and H to reduce the effects of buoyancy of heated, light air at the surface. The algorithm started with finding a sets of “hot” and “cold” pixels. It assumed that LE equals to 0.05 and 1.05 times reference ET (ET_o) for “hot” and “cold “pixels at the beginning of the iteration, and then calculated dT_hot and dT_cold (Equation (4)), followed by calculating a and b, and then updating dT, H, and ${\mathrm{r}}_{\mathrm{ah}}$. We used 5% of reference ET, instead of 0, to initiate the LE for “hot” pixels, considering the possibility of soil moisture content in bare soils shown by [17]. We tested the possible impact of using ET_o for internal calibration, and found that using 1.05 as a constant led to the closest agreement between observed and estimated ET. The Monin–Obuknov length (L) was used to define the stability conditions of the atmosphere during iteration in METRIC [17]. The iteration stops when the value of ${\mathrm{r}}_{\mathrm{ah}}$ remains stable [44] (see Figure S1). After the convergence of the hot and cold pixels, iterations are also needed to converge the surface energy balance components. Dhungel et al. [47] demonstrated the improved behavior of surface energy balance convergence in the low wind speed, which is particularly important to reduce uncertainty in surface energy balance components due to the lack of the convergence. In our study, we set that the iteration process of H calculation stops when the coefficients of a and b remain stable. Our test showed that it typically takes 5 or 6 iterations to reach convergence while low wind conditions required more iterations. A simplified flowchart of METRIC approach is shown in Figure 3.

In METRIC, ET of each pixel at Landsat overpass time (10:30 a.m.) is calculated from Equation (7), and then the fraction of reference ET (ET_oF) is calculated by dividing instantaneous ET by hourly reference ET_o (Equation (8)), obtained from Belridge CIMIS station. It assumes that the value of ET_oF at Landsat overpass time is the same as the average value of 24 h during Landsat overpass day [17,26]. Finally, the daily ET at each pixel is calculated as a product of ET_oF and daily reference ET_o (Equation (9)), which is the accumulative 24 h reference ET_o at the station during the dates when Landsat satellite overpasses every 16 days and when imagery is not obscured by cloud.

$${\mathrm{ET}}_{\mathrm{inst}}=3600*\mathrm{LE}/\left({\mathsf{\lambda }*\mathsf{\beta }}_{\mathrm{w}}\right)$$

(7)

$${\mathrm{ET}}_{\mathrm{o}}\mathrm{F}={\mathrm{ET}}_{\mathrm{inst}}/{\mathrm{ET}}_{\mathrm{O}\_\mathrm{inst}}$$

(8)

$$\mathrm{ET}\_\mathrm{daily}={\mathrm{ET}}_{\mathrm{o}}{\mathrm{F}*\mathrm{ET}}_{\mathrm{O}}\_\mathrm{daily}$$

(9)

where ET_inst is the instantaneous ET in mm/hour, 3600 is used for converting from seconds to hours, λ is the latent heat of vaporization in J/kg, β_w is the density of water, ET_{o_inst} is the hourly reference ET of grass at the time of image, ET_daily is the 24 h actual ET, and ET_o_daily is the accumulative 24 h reference ET_o on the date of image acquisition.

ET-based irrigation scheduling requires continuous daily ET estimate especially during growing season, so that ET losses could be replaced by applying the proper amount of water at the right time to meet the plant water demands. We applied a cubic spline interpolation to fill in the gaps of daily ET_oF using the Landsat-based estimates of ET_oF at the four most adjacent dates when Landsat clear sky observations were available. A cubic spline was first fitted to generate a smoothed ET_oF curve based on the non-linear curve nature of the temporal dynamics of fraction of reference ET [44]. The complete time series of daily ET_oF was then derived from the splined curve, and used in combination with the continuous daily reference ET_o at the Belridge CIMIS station to derive daily ET estimate for every single day. Monthly, seasonal, and annual ET estimates were finally aggregated from daily values. The temporal interpolation was only done during the growing season from April to September, since there were too many missing data during October to March winter rainy season.

3.3. Automated “Cold” and “Hot” Pixel Selection

In order to estimate sensible heat flux (H) accurately, two “anchor” pixels, named “cold” pixel and “hot” pixel must be selected to determine the relationship between Ts and dT in the METRIC approach. “Cold” pixels represent well-watered, full-cover vegetation areas within an image, and thus have relatively higher NDVI and relatively low Ts. The “hot” pixels, namely dry and bare soil, on the other hand, usually have low NDVI but high Ts. The selection of these anchor pixels typically involves some training and manual image interpretation, and is often arbitrary. To make the process more consistent, we developed a method to select those pixels automatically based on the statistics of NDVI and Ts maps, following and integrating similar concepts proposed by previous studies [48,49,50,51]. The cumulative histogram of NDVI and Ts were first generated. The values of NDVI at the point of 5% (95%) and the values of Ts at the point of 95% (5%) in the accumulated percentage curves were used in combination as thresholds for identifying hot (cold) pixels. To allow for flexibility, we extracted all pixels falling both within ±0.01 of the NDVI thresholds and within ±0.5 °C of Ts threshold, and averaged over these two sets of pixels to calculate the mean Ts for “hot” and “cold” pixels (Figure 4). We tested imageries in different seasons and found that the resulting T_s − dT relationships with this automatic selection method were similar to those from manual selections.

3.4. Evaluation of Evapotranspiration Estimate

To assess the uncertainties of METRIC approach in estimating ET over the almond orchard, we compared the METRIC ET with field measurements. Comparisons were made for instantaneous, daily and monthly ET. Considering the fetch area of 50 by 50 m of the installed flux tower in the study area, we averaged the METRIC ET over 2 by 3 pixels centered at the flux tower location. Four statistics were reported, which are mean bias, mean relative difference (MRD), root mean squared error (RMSE) and determination coefficient (R²). The equations of bias, MRD and RMSE are shown as follows:

$$\mathrm{Bias}={\mathrm{ET}}_{\mathrm{METRIC}}-{\mathrm{ET}}_{\mathrm{Tower}}$$

(10)

$$\mathrm{MRD}(\%)=\frac{\left|{\mathrm{ET}}_{\mathrm{METRIC}}-{\mathrm{ET}}_{\mathrm{Tower}}\right|}{{\mathrm{ET}}_{\mathrm{Tower}}}*100$$

(11)

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{ET}}_{\mathrm{METRIC}}-{\mathrm{ET}}_{\mathrm{Tower}}\right)}^2}{\mathrm{n}}}$$

(12)

where ET_METRIC is METRIC-based ET estimate at Landsat 5 TM and 7 ETM+ overpass date and ET_Tower is ground-based ET measurement at the same day from micrometeorological tower. We also examined the ability of METRIC approach to capture the inter-annual and spatial variation of ET.

3.5. Simplified Empirical Approaches for Evapotranspiration Estimate

Many growers do not have the capacity to implement sophisticated remote sensing based algorithm for their irrigation management, especially on a daily basis. We thus also aimed to develop easy-to-use and cost-effective approaches for growers to estimate ET with easily accessible data such as NDVI and weather data in this study. We took advantage of the spatially explicit ET estimates derived from Landsat data, and examined the main factors controlling spatial and temporal ET variation using both temporal and spatial data. We quantified the relationship between METRIC ET, as a dependent variable, and NDVI, which potentially can be relatively easy to derive from imagery by drones, and a suite of environmental variables that growers have easy access to. Multivariate statistical models were built using three years of spatial data. The univariate linear and nonlinear relationships between ET and the explanatory variables were first examined, and then various combinations of linear and exponential terms of key ET driving factors, including NDVI, net radiation, and vapor pressure deficit, were explored. The parameters were optimized using “lsqcurvefit” in MATLAB to minimize the difference between predicted and METRIC-estimated ET. The models were validated with the rest year of both spatial ET data and field measured ET as well.

4. Results

4.1. Evaluation with Flux Tower Data

We compared the instantaneous METRC-based ET estimate using Landsat 5 TM and 7 ETM+ data, averaged over a 2 by 3 Landsat pixel window centered at the flux tower location, with ET values during 10:00 a.m. to 11:00 a.m. derived from flux tower measurements during Landsat overpassing dates. Aggregated ET estimates at daily, monthly, and annual time scales were also compared with corresponding field measurements in order to validate the performance of METRIC approach when it was applied to map ET at different time scales over an almond orchard in the southern San Joaquin Valley.

4.1.1. Instantaneous and Daily ET on Landsat Overpassing Dates

The instantaneous ET during Landsat overpassing time around 10:30 a.m., directly estimated through METRIC approach, agreed well with field measurements in general, with a R² of 0.74 and RMSE of 0.11 mm/h, over a total of 46 Landsat overpassing cloud-free days from 2009 to 2012 (Figure 5a and

Table 1. Accuracy of instantaneous and daily METRIC ET estimates in comparison with flux tower measurements during Landsat overpass dates in 2009–2012.

		Instantaneous ET				Daily ET
Year	Number of Days	R2	RMSE (mm/h)	Bias (mm/h)	MRD(%)	R2	RMSE (mm/Day)	Bias (mm/Day)	MRD (%)
		All clear-sky Landsat overpassing dates
2009	14	0.82	0.11	0.00	27.91	0.87	0.93	−0.02	31.82
2010	11	0.45	0.12	0.02	16.70	0.76	0.91	−0.32	11.69
2011	10	0.89	0.12	0.04	28.47	0.98	0.76	0.37	27.31
2012	11	0.90	0.10	−0.06	13.02	0.96	0.49	0.27	13.67
Total	46	0.74	0.11	0.00	21.79	0.87	0.80	0.06	21.68
		Growing season (April–September)
2009	9	0.69	0.05	−0.02	6.38	0.78	0.56	−0.16	7.56
2010	9	0.86	0.09	0.05	12.64	0.89	0.48	−0.06	6.46
2011	8	0.79	0.11	0.01	15.93	0.94	0.54	0.13	8.31
2012	8	0.60	0.11	−0.06	10.80	0.88	0.54	0.26	9.73
Total	34	0.64	0.09	0.00	11.32	0.87	0.53	0.03	7.96

). When aggregated at daily time scale, a better agreement was found, with a R² of 0.87, RMSE of 0.80 mm/day. Higher accuracy was found during the growing season, e.g., R² remains the same, but RMSE was reduced to 0.53 mm/day and mean relative difference (MRD) lowered to 8.0% (Figure 5b and Table 1).

The METRIC approach captured similar seasonal dynamics of daily ET as the field measurements (Figure 6 and Figure S2). For example, estimated ET varied from 1.61 mm on 3 December to 7.72 mm on 28 July in 2009. Similarly, the measured ET was lowest in January and December, increased gradually until June and July with a maximum of 9.15 mm/day on 20 June, and then declined again, reaching a minimum ET of 0.17 mm/day on 29 December. The large day-to-day fluctuation of field measured ET mostly followed that of the reference ET_o. The mean daily METRIC-estimated ET was 5.10 mm/day, averaged over a total of usable 14 Landsat overpassing days in 2009, which was similar to that from the field measurements (5.12 mm/day) and about 7% higher than the mean daily reference ET_o of 4.80 mm/day. A high bias of the METRIC approach was found in winter in both 2009 and 2010.

4.1.2. Fraction of Reference Evapotranspiration (ET_oF)

The daily time series of ET_oF derived from flux tower measurements showed some large day to day fluctuations, such as peaks (ET_oF > 1.30) on 10 April, 1 May, and 14 September in 2009; 21 April and 17 May in 2010; 8 April, 4 June, and 11 September (ET_oF = 0.54) in 2011, and 13 April in 2012 (Figure S3). The data spikes may be caused by irrigation events or by station maintenance. In addition, main spikes occurred during spring when a local high inversion fog and cloudiness is common and might cause temporal variability of the radiation balance. The METRIC-based ET_oF estimates could not capture this high frequency variation, limited by the 16-day revisiting time of the satellites and cloud cover. We therefore calculated a 16-day running average of field-based ET_oF as the ratio of 16-day running average of measured daily ET divided by 16-day running average of daily ET_o (Figure 7). The complete time series of METRIC daily ET_oF was obtained using cubic spline interpolation from the available ET_oF estimates during Landsat overpassing days, and was then used to calculate the 16-day running average of METRIC ET_oF as show in Figure 7.

At a 16-day time scale, the average values of estimated ET_oF were close to those of measured ET_oF, with a MRD of 6.29% and a small bias of −0.01 from April to September in 2009–2012. Field-based ET_oF was mostly low from January to March, increased significantly from the beginning of April (ET_oF = 0.92) to August or the middle of September (ET_oF = 1.30), and then started to decrease at the end of September in general (Figure 7). For example, ET_oF increased from below 0.6 in January to 0.92 in April, reached its maximum at 1.30 in September, and decreased to 0.6 by the end of year in 2009. A similar trend was observed in 2012, when ET_oF was 0.91 in the beginning of April and reached 1.21 in the middle of September. The ET_oF during Landsat overpassing days followed this similar trend.

4.1.3. Continuous Time Series of Daily and Monthly ET

Daily ET was calculated by Equation (9), as a product of daily ET_o and temporally interpolated daily ET_oF, during days with clouds or without Landsat overpass. We here focused on April to September since there were too many imageries obscured by cloud in winter season and thus the interpolation of ET_oF was expected to result in larger uncertainties. The complete daily time series of METRIC-based ET estimates followed closely with those from field measurements (Figure 8). Overall, METRIC daily ET agreed well with measurements, with a R² of 0.85, a RMSE of 0.56 mm/day, a bias of −0.05 mm/day, and a MRD of 8% for a total of 732 days during 2009–2012 growing seasons (

Table 2. Comparison between complete time series of daily ET estimated by METRIC and measured by the tower from April to September (n = 183 days per year) during 2009–2012.

Year	R2	RMSE (mm/Day)	Bias (mm/Day)	MRD (%)
2009	0.84	0.57	−0.28	7.09
2010	0.87	0.58	−0.06	8.27
2011	0.86	0.53	0.12	8.32
2012	0.83	0.56	0.04	8.31
Total	0.85	0.56	−0.05	8.00

). Comparing with the daily ET estimate during Landsat overpass days only (see Section 4.1.1), the overall accuracy was not decreased after temporal interpolation, partly because all images used in this paper were of good quality (cloud-free) and ET_oF had a small range during growing season. Mean daily ET from METRIC was 6.36 ± 1.24 mm/day averaged during April–September (n = 183 days) in 2009, close to the flux tower measured ET (6.63 ± 1.23 mm/day). Both METRIC and field measurements showed a gradual decrease in mean daily growing season ET in 2010 and 2011, e.g., 5.99 ± 1.42 mm/day in 2010 and 5.73 ± 1.20 mm/day in 2011 from METRIC estimates, and then a slight increase to 6.16 ± 1.31 mm/day in 2012. Similarly, mean daily ET from flux tower measurements decreased to 6.05 ± 1.59 mm/day in 2010 and 5.61 ± 1.36 mm/day in 2011 and increased to 6.13 ± 1.33 mm/day in 2012.

When aggregated at monthly time scale, METRIC ET estimates were in good agreement with flux tower measurements during 2009–2012 growing seasons, as shown by a high R² of 0.90 and a relatively small RMD of 9.68% (n = 24, Figure 9a). The RMSE was 12.08 mm/month and there was only a slight bias of −1.40 mm/month. Both METRIC and flux tower measurements showed a peak ET in July: METRIC ET varied from 207.7 mm/month in 2011 to 232.8 mm/month in 2009, and flux tower measured ET varied from 215.4 mm/month to 243.5 mm/month, respectively. METRIC-based ET was much lower in April, ranging from 120.3 mm/month in 2010 to 156.1 mm/month in 2009, similar with flux tower measured ET.

Satellite-based monthly ET estimates captured similar inter-annual variation as observed by flux tower measurements (Figure 9b). Mean annual measured ET varied from 1027 mm/year in 2011 to 1214 mm/year in 2009 summarized over April–September growing seasons. Similarly, METRIC monthly ET also showed the largest ET across all months in 2009, resulting in a mean annual ET of 1163 mm/year, and a 10% decrease in 2011. Monthly ET from June to September was much lower in 2011 than the other three years, as shown by both METRIC and flux estimates. This was consistent with the lower ET_oFs in 2011 (Figure 7). In addition, ET_o in 2011 was also lower than that of the other three years from April to September.

4.2. Spatial Pattern of Evapotranspiration

Spatial distribution of daily ET at the field-scale is essential to understand the water distribution pattern over a whole orchard to irrigate rationally, especially during crop growing season. Maps of spatial daily ET and main driving factors at a 30-m resolution (676 pixels) are shown in Figure 10, for selected dates representing different growth stages of almond trees. ET within the orchard follows similar spatial patterns of ET_oF, NDVI, and net radiation. Roads, either non-vegetated or sparsely vegetated, as indicated by low NDVI, had lower net radiation and available energy than trees, and smaller ET_oF further contributed to lower ET over sparsely vegetated areas, compared with the rest of the orchard.

Throughout the study period, ET of almond trees was relatively uniform, because the main driving factors, i.e., ET_oF and R_n calculated at Landsat overpassing time, had only slight variations within the orchard. Overall, the coefficient of variance (CV = STD/Mean) of METRIC-based daily ET estimates over the whole almond orchard varied between 3.63% and 10.52% and the CV of monthly ET ranged from 2.60% to 6.24% during the study period (Table S2). The CV of daily ET on 29 June 2010, for example, was 5.52%, with a mean ET of 7.39 mm/day. ET_oF varied between 0.77 and 1.13 and had a mean of 1.04 ± 0.06 (CV: 5.52%), while spatial variability of R_n was much smaller, with a CV less than 2% on the same day. The CV of NDVI was around 8% in general.

The difference between almond trees and the road were minimal during non-growing season, i.e., November and December (not shown), when almond trees go through a period of winter dormancy. The CV increased as the plants developed and then decreased again to dormancy. These findings were consistent with Couvreur et al. [37], which showed that the CV of tree-scale ET was 5–7% and increased as the average soil water stress was larger during 2009–2012 in the same study area.

4.3. Factors Controlling Evapotranspiration

To better understand the main factors controlling spatial and temporal ET variation, correlations between ET and its driving factors were studied using both temporal and spatial data. For each individual pixel in the orchard, the correlations of daily ET with NDVI and R_n were computed for Landsat overpass dates (n = 46 days) during 2009–2012. Figure 11 shows that R_n was a dominate driver in determining day-to-day variation in ET, with an average correlation of 0.96, indicating that the majority of daily ET variation among seasons and years can be captured by R_n. Time series of NDVI was less correlated with that of daily ET, but was still highly correlated, with an average correlation of 0.79 and a range of 0.62–0.85 throughout the orchard (Figure 11b).

The spatial variability of daily ET was also mostly controlled by R_n variation across the orchard for any given Landsat overpass day, especially during April–September growing season, as shown by a typical correlation of 0.96 (Figure 12a). The average growing season correlations between daily ET and NDVI were 0.70, 0.60, 0.64 and 0.87 for 2009–2012, respectively, much lower than those with R_n (Figure 12b). The dependence of daily ET on spatial distributions of R_n and NDVI decreased during non-growing season (Figure 12).

4.4. Empirical Method for Estimating Evapotranspiration

R_n and NDVI were both closely correlated with daily ET temporally and spatially as shown in 4.3. We thus explored the feasibility of using NDVI and R_n in combination with daily reference ET_o, vapor pressure deficit (VPD), and average daily wind speed (U) to develop empirical and easy-to-use methods for daily ET estimation. Spatial data from 2009 to 2011 were used to build a suite of univariate and multivariate models and 2012 data were used for validation and uncertainty assessment. The ET results from the empirical methods were also compared with ET tower measurements. Equation (13) was chosen after model building and testing with all available spatial and temporal data:

$$\mathrm{ET}=\mathrm{ETo}*\left(0.3108*\mathrm{NDVI}–0.0059*\ln \left(\mathrm{VPD}\right)+0.8315\right)$$

(13)

This statistical-based equation indicates that ET_o and NDVI were positively related with ET as expected. VPD, however, negatively regulates evapotranspiration here, possibly as a proxy of water stress on ET_oF [7]. It predicted well the temporal dynamics of evapotranspiration at the tower site during 2012, with a R² of 0.94 and a RMSE of 0.55 mm/day (n = 11 days) when comparing with flux tower measurements, and with a R² and RMSE of 0.98 and 0.27 mm/day (n = 7436 pixels) when compared with METRIC-based ET estimates during 2012 (Figure 13).

We also found that reference ET_o had an exponential relationship with R_n (R² = 0.95) and therefore reconstructed Equation (13) and optimized the parameters to obtain the second empirical algorithm as shown below:

$$\mathrm{ET}=\left[0.9405*\exp \left(0.0027*\mathrm{Rn}\right)\left]*\right[0.1816*\mathrm{NDVI}+0.2672*\ln \left(\mathrm{VPD}\right)+0.8913\right]$$

(14)

The ET estimates from Equation (14) were highly correlated with both tower measured ET (R² = 0.94, RMSE = 0.71 mm/day) and with 2012 METRIC ET estimates (R² = 0.99, RMSE = 0.28 mm/day) (Figure 13). This empirical method is of particular value to estimate ET, when ET_o is not available, since net radiation can be estimated with satellite remote sensing observations. Similar to Penman Monteith concept, VPD here represents more the drying power of an atmosphere and thus atmospheric demand for water.

5. Conclusions and Discussion

The almond acreage in California has grown rapidly in the past 10 years, raising concerns about competition for the limited water resources in California. Information on crop ET at the field-scale is crucial for irrigation management at the farm level and for water planning at the watershed and state level. A robust and cost-effective approach for crop ET mapping is needed to improve sustainable water management. METRIC, a satellite-based surface energy balance method, has been widely used to map crop evapotranspiration (ET) of field crops, but less is known about how it performs for orchards.

In this paper, we assessed the accuracy of the METRIC approach, driven mainly by Landsat satellite data, in estimating daily and monthly ET over an almond orchard in southern San Joaquin Valley, California during 2009–2012 when field measured ET data were available. METRIC-based ET estimates agreed well with flux tower measured ET both instantaneously and at daily time step. High accuracy was found during the growing season (R² = 0.87 and RMSE = 0.53 mm/day). Monthly ET estimates captured similar seasonal ET pattern, increasing from April and peaking in July, and were in good agreement with field measurements with a R² of 0.90 and a RMSE of 13.69 mm/month during 2009–2012 growing seasons. ET estimated by METRIC varied from 119.56 to 226.46 mm/month averaged over three years, and measured monthly ET varied from 121.06 to 243.54 mm/month. Both METRIC and flux tower measurements showed the highest ET in 2009.

Continuous daily fraction of reference ET (ET_oF) was derived from the METRIC-estimated ET_oF during available clear sky Landsat overpass dates using a cubic spline interpolation technique. Overall, the estimated 16-day ET_oFs were close to those derived from field measurements, with a RMSE of 0.08 and a MRD of 6.29%. The interpolated ET_oF curve, however, was unable to capture the high frequency variation due to irrigation events and other factors as revealed by the flux tower measured ET.

We also examined the dominant drivers for temporal and spatial variation in crop water use. METRIC daily ET estimates were highly correlated with R_n and NDVI, both temporally and spatially. The correlation was 0.96 on average between time series of ET and R_n and 0.79 between temporal ET and NDVI (n = 46 days). The correlations between spatial ET and those of R_n and NDVI were 0.96 and 0.70 during 2009–2012 growing seasons, respectively. We further developed two empirical methods to estimate daily ET, using NDVI, R_n, ET_o and VPD as key input data. The method driven by NDVI, R_n, and VPD showed a slightly lower accuracy in mapping ET, with a RMSE of 0.71 mm/day compared with the 2012 flux tower measurements, and a RMSE of 0.28 mm/day compared with the METRIC ET estimate. The empirical method driven by NDVI, ET_o, and VPD showed a lower RMSE (0.55 mm/day) than that of Rn-based method, when compared with the 2012 flux tower measurements.

Overall, the METRIC approach driven by Landsat observations and CIMIS station data can estimate evapotranspiration (ET) in almond orchards with reasonable accuracy during growing season, and also identify both temporal and spatial variability in a cost-effective way. This is consistent with previous assessment of METRIC ET estimates in other crop types. Mean absolute differences between ET estimates from METRIC and field measurements ranged from 16% over irrigated meadow in the monthly scale [28] to 20% in olive orchards in the daily scale [35]. The error was reduced to 4% when aggregated to the seasonal scale [28]. The RMSE varied from 1.1 mm/day in spring wheat fields [32] to 1.9 mm/day over cotton [36].

A continuous daily ET mapping is critical for growers in their efforts to apply the right amount of water at the right time. ET estimate is limited by the 16-day repeat cycle when relying only on Landsat satellite observations, and uncertainties from mathematical interpolations can increase when there is no clear sky imagery during entire month due to cloudy and rainy weather, e.g., during winter and early spring in California. An important next step would be to combine and fuse observations from other sensors such as Sentinel-2 and VIIRS to increase the temporal frequency of remote sensing based ET estimate. The Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM), for example, showed the success and promises of using data fusion to derive Landsat-like daily ET estimates at the field scale [10,52]. Other methods have also been developed to improve the accuracy of three hourly or daily ET mapping, including a combination of the feedback model (GG model) with SEBAL [53]. As discussed earlier, we utilized ET_o to be congruent to the instantaneous and daily ET calculation, in order to interpolate ET_oF between the overpass dates using the spline method. Earlier studies indicate that ET_rF needs to be adjusted in particular crops like olive orchards which behave similarly to the grass-based coefficients [34,35] because of the crop stress, physiological and micrometerological variations. To address this problem, Dhungel et al. [54] developed a backward-averaged iterative two-source surface temperature and energy balance solution algorithm (BAITSSS), incorporating the crop and soil stress components using soil water balance, to capture the impacts of the irrigation and precipitation between the satellite overpass. They demonstrated improved robustness compared with mathematical and regression based interpolations.

Our study here showed that the METRIC ET estimate using well-watered pasture as a reference crop worked reasonably well, although it may not be the case for other crop types. It is always necessary to examine the impact of selection of anchor pixels and the reference ET on actual ET estimate of any untested cropping systems. Further studies are also needed to quantify the ET field measurement errors and investigate the METRIC performance in younger orchards and in more homogeneous orchards under the effects of various environment stressors such as increasing salinity in soil and deficit irrigation.