Remote sensing of evapotranspiration for irrigated crops at Yuma, Arizona, USA

A satellite-based vegetation index model that tracks daily crop growth and evapotranspiration (ET c ) is developed, tested, and validated over irrigated farms in Yuma irrigation districts of Arizona and California. Model inputs are remotely sensed normalized difference vegetation index (NDVI) images, crop type maps, and local weather. The utility and novelty of the model is a more accurate assessment of ET c than currently provided by the US Bureau of Reclamation ’ s evapotranspiration modeling system. The model analyzes NDVI time series data from the European Space Agency ’ s Sentinel-2 satellites using the Google Earth Engine, constructs FAO-56 style crop growth stages from NDVI, and then estimates daily ET c using pre-defined crop coefficients (Kc) and grass reference evapotranspiration (ET os ). Four crops were selected to test and evaluate model performance: short-season broccoli, mid-season cotton and wheat, and perennial alfalfa. Comparison of model results showed that Recla-mation reports overestimate alfalfa and wheat ET c by 21 – 25%, cotton ET c by 6%, and underestimate broccoli ET c by 21%. Variability resolved by the model ranged 6 – 18% of median total ET c . Comparison of model results with those obtained from 13 eddy covariance sites showed validation discrepancies ranging 1 – 14%: average total actual ET c differences were 12, (cid:0) 14, 78, and 87 mm/season, respectively, for alfalfa, broccoli, cotton, and wheat. The wide availability of Sentinel-2 data, collected every 5 days or less, and the rapid processing via Google Earth Engine make the vegetation index model implementation fast and practical. Its accuracy and ability to resolve ET c for every field would benefit the Reclamation water accounting system and provide valuable consumptive water use data for any Colorado River stakeholder.


Introduction
The Lower Colorado River is a critical resource for irrigated agriculture where consumptive use is approximately 3.7 km 3 /year (~3 M ac-ft/year) of water, applied over 335,890 ha (830,000 ac) (Reclamation, 2014).Its availability makes possible a regional economy worth billions of dollars (James et al., 2014).The U.S. Department of Interior, through the U.S. Bureau of Reclamation (Reclamation), manages river resources, and by U.S. Supreme Court Decree (547 U.S. 150 (1963, consolidated 2006)) is required to account for annual water use.The purpose of this study is to introduce an approach that could help improve these water accounting obligations.To better understand why, some background on official reports and techniques is presented.
Reclamation primarily collects, summarizes, and reports water use data in two documents.In one, the 'Colorado River Accounting and Water Use Report' (Water Report, available from https://usbr.gov/lc/region/g4000/wtracct.html), water users in Arizona, California, and Nevada are tabulated with their corresponding uses: water diversions, measured returns, unmeasured returns, and consumptive use The Evapotranspiration Report provides crop evapotranspiration (ET c ) by utilizing crop coefficient curves defined by crop coefficients (K c ) adapted from FAO-56 and ASCE (Allen et al., 1998;Allen et al., 2005), and by crop growth stages.While two K c versions are described in FAO-56, only the single-coefficient model is used, as shown in Eq. 1: where ET os is grass-referenced Penman-Monteith evapotranspiration as described in Allen et al. (2005).Use of ETos, instead of ETo, is done to avoid confusion and ambiguity with Arizona and California weather data sets, where ETo sometimes does not represent FAO-56 grass referenced evapotranspiration (see https://extension.arizona.edu/sites/extension.arizona.edu/files/pubs/az1324.pdf for details).The K c values vary with crop type and are constructed from three reference K c values representing the initial bare soil condition, the mid-season full-cover condition, and end-of-season.Actual crop growth, as represented by fraction cover, usually varies smoothly from day-to-day, but the FAO-56 simplifies growth patterns by partitioning a season into four stages: initial (L INI), development (L DEV), mid-season (L MID), and end-of-season (L END).K c values are constant within the L INI and L MID stages and are linearly connected through the L DEV stage.The final growth stage coefficients are defined by linearly connecting the K c L MID reference value to the K c L END reference value (Fig. 1).Selection of the appropriate reference K c values for the Evapotranspiration Report is determined from crop class maps obtained from remote sensing image data, primarily from Landsat (www.usgs.gov/core-science-systems/nli/landsat)and airborne photogrammetric surveys (USDA National Agricultural Imaging Program [NAIP]).K c values are provided for all major crops and are updated from time-to-time (e.g., Wright, 1982, Mhawej et al., 2021, Pereira et al., 2018, Pereira et al., 2021).
Growth stages are less commonly updated.While there is extensive literature describing how crop growth can be monitored with remote sensing, particularly with the Normalized Difference Vegetation Index (NDVI, Fischer, 1994), the reports focus on crop phenological monitoring, and do not specifically define all four FAO-56 stages.Illustrative examples are provided by Seo et al. (2019), who employ logistic models for tracking corn and soybean, and by Yang et al. (2013), who create web-based tools to identify start, end, and peak of growing seasons.
Nevertheless, it is important to find ways to specify appropriate L values because they are sensitive to local environmental conditions and farming practices.Accordingly, use of generalized growth stages will generate inaccurate time-interpolated K c values.This pitfall is emphasized in p 104): "Lengths of crop development stages provided in this table are indicative of general conditions but may vary substantially from region-to-region, with climate and cropping conditions, and with crop variety.The user is strongly encouraged to obtain appropriate local information."Obtaining localized crop development growth length stages, however, is costly and time consuming, and determining stage lengths when considering differences in crop varieties, crop market conditions, and environmental factors is confounding.Field observations, reports, surveys, and local knowledge must be compiled over several years to yield meaningful growth curves.
One possible option to constrain crop development stage lengths is to use growing degree days (GDD).Use of GDD tends to focus on estimation of time of maturity (Wang, 1960).Studies consider all major crops, examples include soybean (Lou et al., 2023), maize (Cross and Zuber, 1972), sorghum (Conley and Wiebold, 2003), cotton (Peng et al., 1989).More recently, GDD have been used with other observations to improve accuracy, such as solar radiation for lettuce (Wurr, Fellows and Suckling, 1988), or use NDVI for crop type mapping (Skakun et al., 2017).The utility of GDD for crop development mapping, however, may be too generalized: its requirement for knowledge of planting dates must be supplied in other ways.Another difficulty is GDD are statistically averaged generalizations of crop development and thus may not represent local conditions.As McMaster and Smika (1988) point out in their wheat study, GDD are sensitive to variations in agronomic factors such as cultivars, row spacing, fertilizer and planting dates.Water stress, unresolved by GDD, alters growth development and may need additional information to adjust development estimates (Mahan et al., 2014).
Recent transformative developments can improve specificity and overcome obstacles to gaining accurate growth stage data.Improved availability, frequency, spectral sampling, and spatial resolution of remote sensing image data now makes possible the observations of vegetation cover at sub-field scales at nearly daily time steps.With the launch of the European Space Agency satellites Sentinel-2A and 2B in the period 2015-2017, 10 m resolution data could be acquired at least every 5 days, meaning that by using NDVI one can distinguish between bare soil, partial cover, and full-cover fields.In its current configuration, the Evapotranspiration Report uses remote sensing data obtained quarterly for crop classification mapping, but with nearly daily NDVI, the approach could be extended to track growth stages without relying on statically defined L values.
A second transformative development has been the creation of cloudbased satellite image processing services, such as Google Earth Engine (earthengine.google.com)and Sentinel-Hub (www.sentinel-hub.com).Previously all image collections, each exceeding 100MB, had to be located on provider servers, downloaded to local computers, and then processed.Cloud-based data processing enables rapid extraction of fieldspecific NDVI data without the burden of downloading the source data.
An additional consideration for improving the Evapotranspiration Report is its potential performance compared with competing evapotranspiration models.In parallel with remote sensing improvements and availability, the accuracy and robustness of evapotranspiration modeling has improved (Fisher et al., 2017).Two main modeling approaches, both distinct from the Evapotranspiration Report methodology, have become established: 1) evapotranspiration is estimated from empirical relationships between vegetation cover and NDVI (e.g., Johnson and Trout, 2012), and 2) evapotranspiration is modeled by combining physical processes using surface energy balance models and land surface temperature data with NDVI (e.g., Norman et al., 1995).Each has advantages and disadvantages, but relative merits have yet to be fully tested.The recently deployed OpenET project (openetdata.org),currently producing monthly ET estimates from the six different models over the Western USA, will help resolve some questions about accuracy, bias, and robustness of the incorporated six different ET models.
In the first modeling approach, NDVI observations are used to empirically connect remotely sensed images with fractional vegetation cover, and then estimate actual ET c (ET c_act ).Examples include SIMS (Satellite Irrigation Management Support; Melton et al., 2012) and VISW (Vegetation Index ET for the US Southwest; French et al., 2018).These approaches use K c values, crop-specific parameters, weather data, and reference evapotranspiration (Monteith and Unsworth, 1990) to obtain periodic ET c .NDVI-based methodologies use freely available 10-30 m resolution, multispectral satellite image data, enabling the generation of daily ET c at field scales.
In the second modeling approach, ET c is obtained from surface energy balance estimates based on environmental physical conditions and empirical relationships connecting soil, water, vegetation, and atmospheric components.Models include Mapping EvapoTranspiration at high Resolution and with Internalized Calibration (METRIC; Allen et al., 2007), Surface Energy Balance Algorithm for Land (SEBAL; Bastiaanssen et al, 2005), Atmosphere-Land EXchange Inverse (ALEXI) flux with Disaggregation approach (disALEXI; Anderson et al., 2011), and Simplified Surface Energy Balance (SSEBop; Senay et al., 2013).Model physics use land surface temperature observations to constrain energy fluxes, and to estimate plant transpiration and bare soil evaporation.Model empirical relations use NDVI data for fractional cover estimation, and sometimes for partitioning of energy fluxes between plants and soil.Because energy balance models incorporate observations of land surface temperature, they are potentially more robust and more geographically general than the first set of NDVI approaches.Further, energy balance models can readily detect rapid evapotranspiration changes due to water stressed vegetation.However, the land surface temperature sensors collect images less frequently, at coarser spatial resolutions than visible-near infrared sensors and require interpolation over longer times between observations.Our study's objectives are to develop, test, and validate a model that improves and extends the existing Evapotranspiration Report while keeping in mind that multiple alternative ET c estimation approaches exist.The study uses new observations and new tools that could significantly improve the Report in a configuration we refer to as the Vegetation Index Model.The Model updates crop development growth stages for every irrigated field without changing existing K c values.The investigation is reported in three steps.In the first step, methods are developed to extract vegetation cover from NDVI images, to convert the resulting values to growth stages, and to transform these stage values to daily evapotranspiration.Sentinel-2 visible-near infrared image data are accessed and processed using Google Earth Engine tools.In the next step, results from growth stage analyses are presented.Sensitivity test results show the relative importance of growth stage metrics compared to Kc values.Then growth stage delineations are reported for four crops: alfalfa, broccoli, cotton, and wheat.These crops represent approximately 15%, 4%, 5%, and 16% of farm acreage (personal communication, synthesis by C. Sanchez using USDA NASS, USBR, and State of Arizona reports) were selected to evaluate the impact of the approach on sampling short and long-season crops.The evapotranspiration results from the Vegetation Index Model are compared with the Evapotranspiration Report results.In the third step, the Model is validated by comparing modeled growth stages and ET c with eddy covariance evapotranspiration, Evapotranspiration Report-based ET c_act , and independently obtained ET c values derived from the vegetation index-based OpenET model, Satellite Irrigation Management Support (SIMS; Melton et al., 2012, Pereira et al., 2020).

Methods
The evapotranspiration study region evaluates the Yuma-Wellton-Mohawk districts in the Southwestern USA (32.6 • N, 114.0 • W, 30 m above sea level, Fig. 2), which is one of six irrigation regions served by the Lower Colorado River.The districts fall within the hot, arid desert category (BWh) of the Köppen-Geiger climate classification system (Beck et al., 2018).Weather data from the University of Arizona's AZMET (https://cals.arizona.edu/AZMET/02.htm)Yuma Valley site (32.7118• N, 114.7041 • W, elevation 36 m above sea level), and estimates of sky radiation from equations presented in Allen et al. (2005), shows the dominance of clear skies throughout the year (Table 1).This region is ~70,000 ha (173,000 acres) in area and consumes approximately 8.7 × 10 8 m 3 /year of river water (~706,000 ac-ft/year).This constitutes ~30% of Arizona's total consumptive use, 2.99 km 3 /year (~2426,000 ac-ft/year), and ~10% of all use in the US Lower Colorado River, 8.72 km 3 /year (~7073,000 ac-ft/year).
The Vegetation Index Model implemented as described below, is composed of four main components (Table 2): 1) data acquisitions, 2) NDVI processing, 3) crop growth stage identification, and 4) ET c assessment and modeling.Component 1 is outlined in Sections 2.1-2.3.Component 2 is described in Section 2.4.Growth stage estimation (component 3) is described in Section 2.5, and details procedures for selecting day of planting and treatment of single-harvest vs. multipleharvest crops.Assessment steps for component 4 are described in Section 2.6.
Updates are done annually following procedures described in Stehman and Milliken (2007).There are over 14,000 sites across the Yuma-Wellton region.For this study we used only ground-verified sites for the period 2018-2020: ~1100 sites, totaling ~8500 ha, each site no smaller than 2 ha (5 acres).The verified sites represent about 12% of all agricultural land in the Yuma region.To assess improvement to the Evapotranspiration Report approach, we made selections that respectively represent episodic, short as well as long season, and widely distributed crops grown during 2018-2020: alfalfa, broccoli, cotton, and wheat (Table 3).

Table 1
Monthly summary for weather observations obtained for the period 2016-2021 at the Yuma Valley AZMET station.Variables include near surface air temperature (Tair), relative humidity (RH), dewpoint temperature (Tdew), wind speed at 2 m above the ground (U2), precipitation (PPT), grass reference evapotranspiration (ET os ), and the ratio of observed to estimated clear sky solar radiation at typical local satellite overpass time of 11AM (R s /R so ).

The Evapotranspiration Report Crop Coefficient Patterns
The Evapotranspiration Report K c values, day of planting (DOP), and growth duration variables (L), as employed in (Jensen et al., 1990;Jensen, 1998;Jensen, 2003) , are presented for the three single season selected crops (Table 4).K c values, reported verbatim, were developed from FAO-56 and adjusted for local weather (Jensen, 2003).The Report does not use the stress coefficient, K c .Coefficients for alfalfa oscillate between 0.368 for the INI stages and 1.104 for MID stages, while the corresponding length (days) for stages are specified as listed in Table 5.
When combined, crop coefficients and growth stage relationships for most crops exhibit two distinct patterns (green lines, Fig. 3).One pattern (broccoli, cotton, and wheat) follows the single-season trapezoidal FAO-56 trajectory, with low K c values at INI, high K c at MID, and lower K c at END.The second pattern, used for alfalfa, is a quasi-periodic series of trapezoidal FAO-56 trajectories (see forage), where the crop is cut multiple times each year.The Evapotranspiration Report specifies 10 cuttings per year for the Lower Colorado irrigation districts.

Modeling evapotranspiration with weather and crop coefficients
Modeled reference evapotranspiration values were needed to assess effects on ET c from FAO-56 growth duration revisions.We used daily Yuma Valley weather station data.Since 2003, the AZMET network has provided two versions of reference evapotranspiration (ag.arizona.edu/azmet/et2.htm).One version, ET0, is similar, but not identical to the CIMIS standard grass-reference (height 0.08-0.15m; Snyder and Pruitt, 1985); it differs in its estimation of net radiation.The other version uses a short reference crop (height 0.12 m), denoted 'ET os ' (Allen et al., 2005).Following Reclamation practice since 2003, ET os has been used for all crops.Its adoption standardizes reference evapotranspiration values for Arizona and California.
Crop coefficient curves were developed from K c values, as listed in Tables 4 and 5, and mapped to daily time steps, using the Evapotranspiration Report (documented in the Appendix to the 2008 edition) and remotely sensed growth phase data.Thus, K c INI was applied to calendar days between 'DOY INI' and 'DOY INI/DEV'; K c MID was applied to calendar days between 'DOY DEV/MID' and 'DOY MID/END', and K c END was applied to the day at 'DOY END'.K c DEV and K c END values were computed by linear interpolation.K c values prior to start-of-season and following end-of-season were set to zero.Daily ET c were computed using Eq. 1, then summed to obtain season total ET c .

Processing remote sensing data
Mean NDVI time series values for every Reclamation-provided polygon were obtained from Sentinel-2 A/B satellite data (Sentinel-2; https://sentinel.esa.int/web/sentinel/missions/sentinel-2).Sentinel-2 consists of two identical satellites in opposed near-polar sun-synchronous orbits with overpass times ~11 AM local time.Images are collected over 12 visible-near infrared bands at resolutions ranging between 10 and 60 m.Acquisitions for any surface location are collected every 5 days.Due to Sentinel-2's wide swath (290 km) and favorable geographical and climate for the Yuma sites, clear-sky acquisitions were better than nominal, averaging ~3 days.When available, surface reflectance data (L2A) were used.For some early 2018 acquisitions, only top-of-atmosphere data (L1C) were available.
Sentinel-2 NDVI values were obtained using the web-based Google Earth Engine (Gorelick et al., 2017).Similar tools that could be used are Sentinel-Hub (www.sentinel-hub.com) and the R package sen2r (sen2r.ranghetti.info).Retrieval of NDVI data from Google Earth Engine was done using Javascript instructions that linked Reclamation-provided crop maps with Sentinel-2 data for the desired observation periods.The instructions selected, filtered, consolidated, and renamed the desired bands (B04 and B08).The cloud mask band, QA60, was set to zero.Output from Google Earth Engine, provided as comma-separated data files on Google Drive folders, were downloaded and analyzed using R-scripts (R Core Team, 2020;RStudio Team, 2020).NDVI outliers were removed and replaced using a median filter (Brock, 1986).The resulting time series were linearly interpolated to daily time steps, then smoothed with a 7-day window convolution filter (using the stats::filter function in R) to reduce artifacts from clouds and atmospheric corrections.The choice of window width was subjective and based on a value midway between overpass intervals nominally spaced 5-10 days apart.Shorter duration windows would improve resolution of NDVI transitions when data collections had high quality but would also increase false transition selections when the collections had poorer quality.Longer duration windows, on the other hand, reduce false selections, but decrease growth stage transition resolution.Acquisition artifacts were mostly minor at Yuma, however, for cloudier climates, more processing would be needed to build clean NDVI traces (Yang et al., 2022).Results were then written to data tables containing field sites, crop type, observation dates, and field mean NDVI values at daily time steps.
It is important to emphasize that there are substantial benefits in time, cost, and required resources in using cloud-processing.For this study, thousands of images, collectively ~88 GB, and tens of thousands of plot polygons had to be accessed for the 3-year water use study at Yuma.Using Google Earth Engine, currently available for use at no cost, tasks were accomplished in minutes while local data storage was under 1 GB.The alternative approach-local processing-requires hundreds of hours of downloading and processing time and hundreds of GB of personal data storage.Time required increases rapidly with slow data connections.Local processing times are exacerbated when needed data must be retrieved from archives; recent changes at dataspace.copernicus.eu,provider for L2A Sentinel-2 data, now place all but the most recent acquisitions into archive storage.

Identifying crop growth phases with NDVI
Daily field scale NDVI data enable modeling of FAO-56 K c curves by mimicking vegetation emergence, development, maturation and harvest stages.Quantifying fractional cover from NDVI is not part of the proposed Vegetation Index Model.The NDVI time series can be partitioned into growth curves for each field either geometrically or quantitatively.
The geometric approaches transform NDVI by finding curve inflections deduced from time series filters, trends, or curves (ΔNDVI/ ΔTime): transitions between INI, DEV, MID, and END can be located in time by short duration changes in ΔNDVI/ΔTime (e.g., the Moving Average Convergence Divergence (MACD), Gao et al., 2020).Variations can be curve-or trend-based (Gao and Zhang, 2021).The advantage for geometric analyses is their greater objectivity; stage partitioning depends on observations and not on set thresholds.However, it can be difficult to distinguish between short-in-time slope changes that represent actual vegetation growth changes and those that represent acquisition or processing artifacts.Due to this limitation, we used the quantitative approach.The quantitative approach analyzes NDVI time series via thresholds (e.g., Jonsson and Eklundh, 2002, Huang et al., 2019, Le Page et al., 2023); phenological transitions are identified by relative changes in NDVI without consideration of curve slopes.The implementation used in this study, denoted by 'QUANT' for convenience, classifies stages for every field by using the local range (min/max) of NDVI and the selection of transition quantiles that separate the four growth stages.The approach is similar to those presented by Jonsson and Eklundh (2002), White et al. (2009), Cong et al. (2012), Ren et al. (2019) and Huang et al. (2019) for start and end of growing seasons, but differs from them by not employing splines, curve-fits, or seasonal periodicities.We chose not to use those options because they tend to overly smooth stage transitions.NDVI partitioning quantifies the duration of the four FAO-56 growth stages and relies upon reference L INI values (specified in the Evapotranspiration Report) to constrain day of planting estimation.QUANT partitions the observed NDVI time-series with quantiles derived from field-specific early-season minima and mid-season maxima (Table 6 and Fig. 4).Minimum and maximum NDVI values are plotted as open circle symbols.NDVI quantiles, corresponding to probability levels of 10%, 90%, and 50% are computed relative to these extremes.Selection of appropriate levels in other studies has been guided by optimizing correlation between thresholds and historical day of planting (Huang et al., 2019), but day of planting values are not generally known, while intervals between day of planting and plant emergence depend upon crop types, environment, and management practices.These intervals are commonly 3-4 weeks for crops grown at Yuma.For the Vegetation Index Model, threshold selection was focused on the identification of the   Jonsson and Eklundh (2002).This correspondence is based on field practice-day of planting is close in time to initial irrigation events-and the fact that soil wetting usually decreases NDVI values.Nevertheless, this brief-in-time NDVI depression is not assured because bare-soil NDVI is already close to its minimum.Thus, a selection constraint is needed.For the Vegetation Index Model, day of planting was chosen assuming it corresponds with the early season NDVI minimum as long it lies within 10 days of the nominal FAO-56 duration.Because the FAO-56 duration values are estimates, use of a time window achieves two objectives: first, it allows NDVI observations to govern day of planting selection as long as those observations are within expectations, and second, it provides a reasonable alternative selection when NDVI-minima derived estimates lie outside expectations.Use of 10 days to govern the expected time interval was not arbitrary but was obtained by analysis of observed and estimated day of planting values for validation sites (presented below, Section 3.3).
Thus, the day of planting determinations for this study were done as follows.First, the day of year corresponding to the early season NDVI minimum is tentatively considered the day of planting.Next, nominal durations of the initial growth stages for the considered crops were determined from FAO-56 tables.These were 35, 50, and 20 days for broccoli, cotton, and wheat, respectively.Alternative days of planting estimates were then determined by subtracting the crop-specific FAO-56 duration from the NDVI-derived day corresponding to the INI/DEV transition.Expected time windows were computed by subtracting and adding 10-day buffers to these alternative days.Last, the tentative days of planting were compared with the time windows and a day of planting selection is made.If the tentative NDVI-based day lies within the window, its value is used, but if it lies outside, the alternative day is used.
The last stage, END, is identified using a 50% quantile-a value chosen to signal end of season while avoiding threshold detection failures that would occur if a lower quantile were used.Note, however, that the end of season for vegetable crops corresponds to final harvest of green plants and not to senescence.For non-vegetable crops such as wheat and cotton, the senescence stage does exist, but it is usually brief because of cropping system practice.At Yuma, precedence is given to higher value vegetable crops and field preparation for them, which means that season lengths for spring and summer crops are commonly terminated by early August.
The steps just described are applicable for vegetable crops, but modifications to the thresholds, and procedures to identify them, are Fig. 4. Example of quantile selection of growth phases for a wheat field in Yuma, 2020, using filtered daily NDVI (heavy blue line).The minimum and maximum values (circle symbols) determine the range for computing quantiles.Horizontal gray lines denote 0-100% quantiles at 10% intervals.The intersection of 10% and 90% NDVI quantiles (bold horizontal red and green lines at NDVI values of 0.19 and 0.80 respectively) with daily NDVI locates growth stage transitions between INI, DEV, and MID.The 50% quantile (yellow line) intersection, during vegetation decline, occurs at the right-most full circle and denotes the last day of the growing season for this field.needed for alfalfa (Table 7).Unlike the single-season crops, alfalfa exhibits multiple NDVI minima and maxima, and has abbreviated INI and END stages.While FAO-56 documentation lists L values for first and subsequent cuttings, use of remotely sensed NDVI can be much more accurate and specific: absent cloudy sky intervals, it tracks all cuttings.However, due to short durations between growth stages-cutting transitions are brief -it is sometimes difficult to resolve them within 5-day remote sensing periodicities.To accommodate this obstacle, the ENDstage duration is set to span the interval between local NDVI maxima and the subsequent NDVI minima.This usually results in END durations of 1-2 days.Local NDVI minima demarcate the start of an INI-stage.
Transitions between INI, DEV, and MID periods are selected as before using 10% and 90% quantiles.Because there are multiple stage transitions within a year, quantile selections are done separately for each cutting cycle (Fig. 5).Minima (circles) and maxima (triangles) are identified within each interval and used to compute quantiles.The alfalfa cutting stage separations are accomplished by using a seasonal trendline that separates local NDVI minima and maxima (red line, Fig. 5).By interpolation of the trendline, the days when the NDVI of the seasonal trendline equals the filtered daily NDVI line can be found.These days demarcate the cutting and full cover time intervals and make possible the identification of NDVI extremes for each cutting cycle.Every field has its own cutting pattern and unique seasonal trendlines must be constructed for each.The approach adopted was to create lowpass filtered NDVI trendlines based on an experimentally determined smoothing window size.Ten NDVI traces over randomly selected alfalfa fields were selected, and trendlines were created using window intervals ranging between 10 and 150 days.Longer time intervals work better for mid-season discrimination, while shorter intervals avoid excessively long taper-on and taper-off periods.Intervals ranging between 50 and 90 days, with constant end-tapers, were suitable for the Yuma samples and the mid-value (71) was selected for all sites.Other sites may require a different filter window width to develop seasonal trend curves.

Eddy covariance data
Eddy covariance data were used to validate the proposed Vegetation Index Model (Table 8).Except for alfalfa sites, day of planting for each site is known-they occurred within one day of eddy covariance deployments.Evapotranspiration values were assessed from 13 sites over 4 crop types: alfalfa, broccoli, cotton, and wheat.Observations were analyzed using EddyPro 7 software in express mode (Fratini and Mauder, 2014;LI-COR Biosciences, 2021).The following corrections were applied: density adjustment (Webb at al, 1980), tilt adjustment with double rotation (Wilczak et al., 2001), block averaging to 30 min, time lag compensation, spike removal, and spectral corrections at low and high frequencies (Moncrieff et al., 2004;Moncrieff et al., 1997).Flux footprint analyses were computed using the method of Kljun et al. (2004).Gaps in weather variables (air temperature, humidity, wind speed, wind direction) were filled using interpolated hourly AZMET data.Missing net radiation values were filled using estimates derived from linear modeling net radiation against AZMET incoming solar radiation.Missing soil heat flux values were filled by replicating values from preceding corresponding hours.Missing sensible heat flux and latent heat flux values were estimated using the online version of REd-dyProc (https://bgc.iwww.mpg.de/5622399/REddyProc;Wutzler et al., 2018).Friction velocity (U*) thresholds were selected using the moving point test (Papale et al., 2006).This tool could not be used for two broccoli data sets because their duration times were < 90 days; for these sites, gaps were filled according to the gap duration.For gaps less than 6 h, H and LE fluxes were linearly interpolated.Longer gaps were filled using local-in-time net radiation and Bowen ratio values.Corrections for changes in heat storage above heat flux plates and photosynthesis were applied (Meyers and Hollinger, 2004).Following Leuning et al. (2012), closure was assessed at 24-hour intervals and enforced using the average daytime Bowen ratio.Quality of eddy covariance results was quantified by regressing daily radiative energy (Rn-G) against advective energy  (H+LE) and reporting coefficient of determination (R2) and regression slopes (b1).Using corrected eddy covariance observations of ET c_act , validity of total season ETc from the Evapotranspiration Report, Vegetation Index, and SIMS models was assessed in three ways.First cumulative ET c obtained from each model was compared with EC ET c_act site-by-site, then summarized by the four crop types.This comparison quantifies absolute and relative ET c estimation bias.Second, cumulative observed ET c values are regressed against modeled values.These regressions summarize model performance with R 2 and root-mean-square (RMSE) statistics.Third, example daily time-series ET c plots illustrate the time and magnitude of model discrepancies.These discrepancies can be localized in time, as might occur during pre-emergent periods when vegetation index models are expected to be less accurate or globalized and thus an indicator of crop coefficient bias error.

Results
This section presents results in four parts.First, the sensitivity of modeled ET c to growth stage delineation is reported for single-harvest crops broccoli, cotton, and wheat.Next, the observed NDVI and growth patterns are reported for alfalfa, broccoli, cotton, and wheat for all Reclamation-verified fields at Yuma.In the third part, a summary of modeled ET c results for the 856 Reclamation-verified sites is presented.Last, ET c outcomes from the model approaches-provided by the Evapotranspiration Report, the Vegetation Index Model, and the OpenET SIMS model-are validated against eddy covariance measurements at 13 sites.

Sensitivity of Modeled ET c to Growth Stage Delineation for Broccoli, Cotton, and Wheat
Using the Evapotranspiration Report K c curves and 2019 AZMET ETos values from Yuma Valley station, the sensitivity of total modeled ET c to changes in L and K c was evaluated for broccoli, cotton, and wheat.Alfalfa results are not presented because partitioning effects are less important than the effects from under-or over-counting cuttings.Inclusion or exclusion of a single cutting can change ET c by more than 150 mm, while changes in L INI, DEV, and MID cause smaller differences, 10-40 mm.
Four L sensitivity tests were conducted, all retaining the nominal length of season: one using the coefficients and growth stages from the Evapotranspiration report as-is but shifted forward or backward in time with respect to its nominal position (DOYshift), and three that varied the durations for the INI, MID, and END stages.Results are presented in the top row of Fig. 6.Effects from all four sensitivity tests were non-linear and seasonally dependent.Changing the start of the growing season (DOYshift) had the most significant effect on changes in estimated total ET for broccoli and wheat, while changes in L INI were most significant for cotton.Twenty-day shifts affected estimated ET for broccoli by 10% (25 mm), for cotton ET by 5% (60 mm), and for wheat by 17% (100 mm).For summer cotton crops, changes in L INI can be relatively more important than DOYshift, which showed that a 20-day deviation from nominal (51 days) changes modeled ET c by 6% (70 mm).
By comparison, sensitivity tests for changes in crop coefficients (Fig. 6, bottom row) showed that accurate K c MID were most important for all three crops, while accuracy for K c INI and K c DEV were much less important for assessing total ET c .A 10% change in K c MID effected a 20 mm (7%) change for broccoli.Sensitivity to changes in recommended K c values depends on the crop.For broccoli, Pereira et al. (2018) proposes a 5% increase in K c MID (1.05-1.10),indicating minimal changes to ET c .For cotton, similar changes to K c MID for cotton would be more important, where seasonal ET c increases by 60 mm.For wheat a 5% increase in K c MID results in a 30 mm ET c increase.
Overall results show that accurate delineation of stages is just as important as accurate crop coefficients.

Crop NDVI and Growth Stage Patterns at Reclamation Verified Sites
Using Sentinel-2 data from 2020 illustrates the NDVI patterns for example crops (gray lines, Fig. 7), and how they correspond to nominal Evapotranspiration Report curves (dashed black lines).The NDVI traces document the range, variability, and timespan for each crop.For the annual crops there is a clear relationship between NDVI and K c .
Alfalfa NDVI patterns (Fig. 7, lower right) are complex, where cuttings occur anytime and do not follow the idealized Evapotranspiration Report pattern (dashed black line).Maximal NDVI is ~0.8 for most times, except for July-August when it drops to 0.6-0.7.This depression likely reflects management practice to reduce summer-time irrigation amounts.For some sites the alfalfa crop ends mid-year, as indicated by horizontal traces where NDVI values are close to zero.
Broccoli NDVI patterns (Fig. 7, upper left), show a wide range of growing seasons with plantings initiated from mid-September (Day of Year 260) through late November (Day of Year 330).The Evapotranspiration Report K c curve represents nominal October planting and falls within the mid-range of the NDVI traces.Full cover NDVI patterns for broccoli, ~0.7, have variable time shapes, with some traces nearly flat while others are rounded.This pattern variability makes stage discrimination for broccoli more difficult than for other crops.
Cotton NDVI (Fig. 7, lower left) documents planting from February to April, which reveals abrupt increases from 0.2 to 0.7 during crop development, while the index is flat-topped during full cover between May and August.The NDVI signals drop in early September, indicative of end-of-season.This pattern is consistent with common Yuma farm practice, where cotton crops must be removed to allow time for prevegetable season leaching fraction irrigations.However, the Evapotranspiration Report K c curve models an atypically longer cotton season.
Wheat NDVI data (Fig. 7, upper right) shows traces almost perfectly horizontal at full cover (NDVI 0.9-0.95),planted during January and harvested in a narrow time interval in May, and closely follows Fall crops (note preceding times with descending NDVI).The corresponding Evapotranspiration Report K c curve suggests an early December planting that is not supported by satellite observations.Differences between QUANT based Vegetation Index Model and Evapotranspiration Report-specified growth stages for the annual crops (broccoli, cotton, wheat) between 2018 and 2020 are shown in Table 9. Listed for each are Evapotranspiration Report and Vegetation Index Model estimates for day of planting and for the four growth stage lengths.The differences in average number of days (Δ) are substantial.Changes in all recorded cases commonly exceed one week, and sometimes one month.The Vegetation Index Model revises day of planting times and, implicitly, times of initial irrigation events, by showing that the actual start of season can be earlier or later than specified by the Evapotranspiration Report: broccoli plantings are 54 days earlier while wheat plantings are later by 34 days.Vegetation Index Model estimates of L DEV are reduced for all three crops, with a major reduction of 43  days for cotton.L MID lengths were increased by 7-12 days for broccoli and cotton and decreased by 24 days for wheat.Considering all growth stages, the Vegetation Index Model revises total season lengths, indicating that Evapotranspiration Report estimates are too long by ~6 weeks for cotton and wheat, and too short by 18 days for broccoli.When combined with sensitivity test results (Fig. 6), the estimates provided by the Vegetation Index Model indicate that significant updates in modeled ET c are needed.Considering growth stage for broccoli (Table 9), day of planting values generated with the Vegetation Index Model were earlier and total L values were greater than the Evapotranspiration Report.For cotton, Vegetation Index Model revisions were in the opposite direction, where days of planting were later, while L totals were less, than provided by the Evapotranspiration Report.Similarly, Vegetation Index Model estimates for wheat sites showed days of planting about one month later, and L total values less than, Evapotranspiration Report estimates.Growth stage results from Vegetation Index Model analyses over alfalfa sites (Table 10) were significant: while the number of days between cuttings were identical for both models, the average number of cuttings per year were 8 for the Vegetation Index Model, in contrast to the Evapotranspiration Report cuttings number of 9-10.
Assessment of remote sensing growth stage pattern retrieval with the Vegetation Index Model for alfalfa (477 sites, 2018-2020) was also done, but using different metrics from the other crops because its stages are short and asynchronous.Instead of tracking transitions between INI, DEV, MID, and END for each cutting event, the durations between cuttings and total number of cuttings per year were tracked.
The first metric considered, duration, is important to estimate because the assumption of a constant value-as is done with the Evapotranspiration Report-can lead to inaccurately estimated seasonal ET c .This is reinforced by the fact that crop management objectives, whether nutrient or yield focused (Orloff and Putnam, 2004), drive the farm-chosen durations, and these may change from year-to-year.Data from 2019 (Fig. 8) illustrates the variability in durations and the general trend where they are longer in the winter than the summer.Shown are the durations and mean day of year for each cutting, coded by colors.The first cutting of the year (red dots) has a mean duration of ~40 days and grows between early January to mid-February.Subsequent cutting events until summer (orange, light green, medium green, turquoise symbols) progressively shorter durations, with a minimum of ~30 days.Cutting events in the fall have progressively increasing durations (medium blue, dark blue, purple, magenta).The median durations (black line) have a similar pattern to data observed at lysimeters in Phoenix (personal communication, D. Hunsaker).For cut durations ranging between 20 and 80 days at Yuma, ET c /cut in the summer ranges ~150-600 mm and ~50-200 mm in the winter.
Nominal durations from the Evapotranspiration Report for alfalfa are 60 days for the first and last cuts of a year, and then decrease to 30 days for other times (Table 10).Results from the Vegetation Index Model over 257 Yuma sites show that durations were nearly the same for most cutting events, except for the first cut of the year, when the Vegetation Index Model observed an average duration of 42 days.
The second metric considered, the number of cutting events per year, is important because the effect of incorrect accounting on seasonal ET c can be large, about 10% per mis-counted event.The Vegetation Index Model observed an average of 8 cuttings per year, while the Evapotranspiration report provides more than 9 per year.Analyses from all fields also show that alfalfa cutting frequency usually ranged between 5 and 9 per year (Fig. 9).Fallowed fields are represented by the zerocuttings entry.

Vegetation Index Model ET c at Reclamation Verified Sites
Vegetation Index Model results of ET c over Yuma Reclamation sites were compared with those returned by the Evapotranspiration Report.Day-to-day patterns of daily ET c use 2020 data (Fig. 10.Red lines indicate ET os while green lines denote Evapotranspiration Report ET c .Daily Vegetation Index Model ET c (VI, black lines) trace site-mean seasonal patterns.Site-to-site variability (blue shading) represents the interquartile range of Vegetation Index Model estimates of ET c .The plots illustrate when and where the Model estimates deviate from the Evapotranspiration Report.Noteworthy instances occur for early season wheat, late season cotton, and most of the year for alfalfa.Wheat and cotton fields are spring to summer crops and show progressive increases in ET c until senescence and harvest; site-to-site variability for these two crops is low at mid-season.Broccoli, a fall-winter crop, has a nearly flat daily ET c profile, with high site-to-site variability for all but the middle of the growing season.Alfalfa, a full-year crop, has a daily ET c pattern that increases with day length, with periodic decreases throughout the year due to cutting events.
Aggregating ET c by crop, year, and model, showed that differences between the Vegetation Index Model and the Evapotranspiration Report model were large for broccoli (+0.21), wheat (− 0.25), and alfalfa (− 0.23), but modest for cotton (− 0.06) (Table 11).Total Vegetation Index Model ET c estimates for broccoli are greater than Evapotranspiration Report by 51-80 mm per year.Vegetation Index Model estimates for cotton are less than Evapotranspiration Report, ranging 30-109 mm per year.Vegetation Index Model results for wheat show that it estimates lower water use than provided by the Evapotranspiration Reportreductions range between 93 and 286 mm.Still larger modeling differences occur for alfalfa (388 and 479 mm), where the Vegetation Index Model returns an average 1478 mm vs. ~1900 mm for the Evapotranspiration Report.These differences are mainly due to corresponding differences in accounting for cutting events.The Vegetation Index Model usually observed 7-10 cuttings annually, while the Evapotranspiration Report specifies 10 for all locations and times (Fig. 9).

Validation of vegetation index model growth stages and actual crop evapotranspiration at eddy covariance sites
Using NDVI data from Sentinel-2, and ET c _act results from eddy covariance observations at 13 sites (Table 8), allows evaluation of Vegetation Index Model accuracy for tracking crop growth stages and for estimating ET c .
Day of planting values obtained from the Vegetation Index Model and eddy covariance data were compared for the single-season cropsbroccoli, cotton, and wheat.Day of planting accuracy when using the early-season NDVI minimum varied widely (Table 12), where estimates ranged from 19 days too early for broccoli and 30 days too late for cotton.When using FAO-56 INI durations, the estimation errors were less, ranging between 10 days too early for cotton to 9 days too late for cotton and broccoli.These smaller differences when using FAO-56 led to the previously noted use of a ± 10 day window as a constraint for the Vegetation Index Model.
The impact of the adoption of a day of planting model constrained by FAO-56 time windows can be assessed using Reclamation-verified sites where the actual day of planting is not known (Table 13).Listed are the fractions of sites where use of an early-season NDVI minimum falls within the FAO-56 window.These fractions range widely, between 0.35 and 0.82, which means that NDVI minima are informative but are often inaccurate.
For the continuous crop, alfalfa, growth stage accuracy was assessed by comparing modeled vs. observed days of cutting (Table 14).Eddy covariance identifications were based on evaluating when daily ET abruptly decreased, while Vegetation Index Model identifications used outcomes from the partitioning of NDVI (Fig. 5).The Model accurately tracked all cutting events at the three validation sites.Average days-ofdetection derived by the Model were different from eddy covariance observations by 2 days or less.
Validation tests of Vegetation Index Model ET c are summarized in Table 15.Shown are crop types, site names, eddy covariance deployment/removal dates, and total season eddy covariance-measured ET c_act .
The last three columns in the table list modeled total ET c from the Evapotranspiration Report (ER), the Vegetation Index Model (VI), and SIMS.Reporting of SIMS values were obtained from the openetdata.orgonline tool and checked locally following methods described in Melton et al. (2012), Pereira et al. (2020), and "The Satellite Irrigation Management Support (SIMS) System User Documentation", 2018-2021 (available online from openetdata.org/methodologies/).To demonstrate that the eddy covariance-based ET c_act data were reasonable, energy balance closure errors were computed and total ET c_act values were compared with literature-reported values.They were generally consistent with other observations.Average closure at daily time steps was 0.82 or better.For broccoli, eddy covariance Fig. 8. Remotely sensed cutting durations vs. day of year for alfalfa grown in Yuma, 2020.Unique colors represent for every field the duration and mean day of year for a specific cutting interval.Intervals 1-9 are plotted as red, orange, light green, medium green, turquoise, medium blue, dark blue, purple, and magenta.For example, a red dot denotes the duration and mean day of year for the first cutting event of a single field in 2018, and an orange dot denotes the same for the second cutting event.The black line and the number at each vertex denote the median duration and day of year for each cutting event.Experimental data from the lysimeter studies at the U.S. Water Conservation Laboratory, USDA/ARS, Phoenix, AZ are shown as a gray line.A.N. French et al. measurements returned an average of 271 mm, compared to 359 mm (post-transplanting, Lopez-Urrea et al., 2009), 220 mm (lysimeter results, Bryla et al., 2010), and > 203 mm (Johnson and Trout, 2012).For cotton, average eddy covariance-based total ET c_act was 935 mm, compared to 918-1033 mm for a Maricopa, Arizona study (French et al., 2015).For wheat, average eddy covariance-based ET c_act was 590 mm, while 478 mm was reported for Maricopa studies (Hunsaker et al., 2005).Eddy covariance-based ET c_act for three alfalfa sites were 1473 mm (YMIDD21-22b1, 301 days), 1149 mm (YMIDD21-22b,220 Fig. 10.Daily modeled ET c for all Yuma instances for four crops in 2020.Shown for each are mean Vegetation Index Model ET c (VI, black), Evapotranspiration Report ET c (ER, green), and grass reference ET (ET os , red).The envelope of results from all fields (blue shading), represents the site-to-site variability of daily ET c .

Table 11
Total median crop ET c for all ground-validated broccoli, cotton, wheat, and alfalfa sites, 2018-2020, as estimated by Evapotranspiration Report and by the Vegetation Index Model.ET c variability estimated by the Vegetation Index Model, as represented by median absolute deviation (MAD), is listed in parentheses, Δ -difference in total median ET c between Vegetation Index Model and Evaporation Report, Rel.Δ -direction and fraction of change in the recorded Δ.   days), and 1809 mm (YMIDD21-22c, 376 days).These were consistent for the first and third sites: model estimates using the IAS model (Snyder and Bali, 2008), returned 1443 mm, and 1922 mm respectively.Inconsistencies between EC and modeled values at YMIDD21-22b could not be explained.The eddy covariance data were checked and adjusted for closure errors and flux footprint outliers and are considered representative.

Table 14
Comparison of Vegetation Index Model (VI) vs. eddy covariance (EC) observations of alfalfa cutting events.For each cutting, the first value indicates the observed day of year from the Vegetation Index Model and the second indicates the day of year estimate from EC data.Fig. 11.Validation results from comparisons of eddy covariance total ET c_act mm (EC ET) for four crops: alfalfa (red), broccoli (green), cotton (blue), wheat (purple).Presented are three ET models: Evapotranspiration Report (ER), Vegetation Index Model (VI), and SIMS.Linear regression are blue, the one-to-one line is plotted in dotted black.
Next, the models were assessed by aggregating results from all four crops and comparing seasonal ET c with eddy covariance observations (Fig. 11).Shown are regression results for each model and the one-toone line (dashed).All linear models have high accuracy with R 2 0.95-0.96(Table 16).However, the Vegetation Index performs better than the Evapotranspiration Report and the SIMS model at ET c values less than 1000 mm as shown by the small (− 4.8) zero-offset term.
Lastly, the validation using ET c daily time series (Fig. 12) highlight reasons for ET c modeling inaccuracies.One lies with the Evapotranspiration Report's wheat growth stages in comparison with satellite derived stages.When viewing daily ET c_act patterns, the Vegetation Index Model and SIMS accurately track eddy covariance measurements for most of the cropping season, but the Evapotranspiration Report erroneously initiates growth prematurely (day of year 360-25 in the new year) and consequently over-estimates early season ET c .Because all three models under-estimate bare soil evaporation, the apparent superiority of the Evapotranspiration Report for seasonal wheat water use is illusory: estimates of high early season ET c result from prematurely advancing INI, DEV, and MID growth stages.A second reason for model inaccuracy, and one that applies to all three models, is their poor estimation of INI-stage evaporation.The early season errors, significant for broccoli (day of year 250-325) and cotton (day of year 50-110), are visually apparent (Fig. 12) by the large differences between observed ET c_act (eddy covariance, EC, black lines) and modeled ET c via ER, SIMS, and VI (respectively, green, blue, orange lines).

Discussion
This study's objectives were to develop, test, and validate the capability of the Vegetation Index Model to extend and improve the existing Evapotranspiration Report methodology for quantifying ET c .A related goal was to consider ET c results against a competitive model, SIMS.Development, testing, and validation were possible because of the availability of three components: validated crop type maps, frequent high resolution satellite data, and software to process the image data.
The wealth of accurate crop cover type maps made possible the matching of NDVI data to hundreds of field sites.These Reclamationsupplied maps provided field-level crop type verification that were assembled from in-person surveys conducted quarterly and made available for a 3-year period.In contrast, studies validating K c values have relatively few sites.For example, Lopez et al. (2009) Ko et al. (2009), and four fields for two years (Hunsaker et al., 2005).For alfalfa, Hu et al. (2020) had treatments in 2018-2019, Benli et al. (2006) investigated one field for three years., and three sites in 1985 were evaluated by (Hunsaker et al., 2002).
Availability of Sentinel-2 satellite data enabled fine temporal scale detection of crop growth changes.Prior to this mission, lower resolution, multispectral data from sensors such as Landsat were available no more frequently than 14 days, and less often with cloudy skies.Deployment of the Sentinel-2 satellite pair has enabled 3-to 5-day observation intervals, providing a capability with worldwide ramifications (Li and Roy, 2017), and facilitated by favorable Yuma weather.High revisit frequencies were especially helpful for alfalfa due to frequent harvesting cycles.
An outstanding component was the computational processing tool provided by Google Earth Engine.Without it, the comprehensive collection, filtering, and reduction of satellite based NDVI would have been slow.Google Earth Engine enabled hundreds of Sentinel-2 images to be combined with crop maps to produce vegetation index time series tables in a few minutes.The high processing speeds mean that annual analyses would be practical for all Lower Colorado irrigation districts.
A quantitative threshold-based methodology, QUANT, was developed and incorporated into the Vegetation Index Model to potentially assist with the Evapotranspiration Report's inability to account for current and future cropping practices.QUANT has similarities to other phenological tracking approaches (e.g., Jonsson andEklundh, 2002, Cong et al., 2012), but differs from them by its use of cloud-processing and constraints specific to transforming NDVI into the four FAO-56 growth stages.With the Evapotranspiration Report, all fields of the same crop type in a particular irrigation district and year are modeled in the same way-identical growth stages and K c values are used, resulting in unrealistically uniform ET c estimates.With the Vegetation Index Model, every field is modeled individually using NDVI-derived growth stages.The resulting availability of remotely sensed ET c makes it possible to evaluate water use every year for all Lower Colorado River irrigation districts.Using crop classification maps and surface weather data, the Model could be readily extended to other irrigation districts.
Testing of the Vegetation Index Model showed that growth stages and ET c values were substantially different from the Evapotranspiration Report (Jensen, 2003), often exceeding 25%.Additionally, the Vegetation Index Model quantified the very important variability of ET c between fields, values not estimated by the Evapotranspiration Report.Two kinds of tests were performed: one evaluated the effect of growth stage beginnings and durations on ET c over Reclamation-verified sites, and the second established relationships between ET c and growth stage variability by conducting sensitivity analyses.
Using remote sensing, comparison tests between models over 856 Reclamation-verified fields-for the period 2018-2020, for alfalfa, broccoli, cotton, and wheat-showed that changes in estimated total ET c values were considerable (Table 10).Vegetation Index Model-derived ET c values were increased by more than 20% for broccoli and decreased by more than 20% for wheat and alfalfa.Absolute ET c differences for wheat and alfalfa were large, where seasonal differences exceeded 100 mm and 400 mm respectively.Total ET c difference for cotton was less, ~6%, because the effect of the Vegetation Index Model's relatively shorter growing season was offset by a shift of peak growth to days with high ETos.
Sensitivity analysis testing of the Vegetation Index Model for fields planted in 2019 (section: 3.1) supported assessments shown in Figs.7 and 13: large growth stage variability occurs within irrigation districts.Because of its use of static K c curves, the Evapotranspiration Report cannot resolve these ET c variabilities.Adoption of the Vegetation Index Model would reduce ET c estimation bias.For the single-harvest crops, broccoli, cotton, and wheat, the analyses showed which growth stages were most important for ET c mapping.In most instances, accurate day of planting was most important, where plausible revisions of 20 days or more can lead to seasonal ET c changes of 25 mm, 60 mm, and 100 mm for broccoli, cotton, and wheat.L INI and L MID were important for broccoli and cotton, while L MID and L END were important for wheat.Crop coefficient uncertainty can be important too, however, sensitivity analyses for Yuma suggest that updates to these are less influential.Updates to K c INI and K c MID are 5% or less (e.g., vegetable crops, Pereira et al., 2018), meaning that corresponding updates to seasonal ET c values would be ~8 mm, 52 mm, and 34 mm, for broccoli, cotton, and wheat respectively.For multiple-harvested alfalfa, ET c is about 200 mm per cutting interval and constitutes about 10% of seasonal ET c .Thus, the number of such periods is the most important growth metric to track for alfalfa.Duration of each interval is less important: revision of durations by 10 days would change ET c by ~50 mm, ¼ of the effect realized by miscounting of events.Validation of the Vegetation Index Model showed that the approach is accurate for alfalfa, broccoli, cotton, and wheat.Using eddy covariance data collected at 13 Yuma sites showed high estimation accuracy, R 2 = 0.95 and RMSE of 99 mm, supporting adoption of the NDVI-driven growth stage partitioning approach.Total ET c values for Vegetation Index Model-derived alfalfa, broccoli, and cotton differed from eddy covariance measurements by less than 8%.Model estimates for wheat  differed by a greater amount, 14%.OpenET-based SIMS estimates for wheat differed by 17%, which suggests that its additional model complexity (e.g., its use of the density factor, K cd ) does not result in superior ET c estimates.The inability to resolve early-stage evaporation by these models was important for cotton and broccoli, but inconsequential for alfalfa and wheat (Fig. 12), suggesting that bare soil durations for these latter crops are brief.
Eddy covariance data also showed high ET c accuracy for the Evapotranspiration Report approach (Fig. 11, R 2 = 0.95).This positive outcome, however, should be considered alongside observations that show its inability to discriminate ET c variability within irrigation districts.Consolidated displays of K c values for the four crops (Fig. 13), show that the field-to-field variability in growth stages (gray lines) is wide in time, on the order of 50 days.The 13 validation sites (red lines) sampled those differences for alfalfa, broccoli, and wheat, and less satisfactorily for cotton.With an exception for broccoli, Evapotranspiration Report K c curves (heavy black lines) are biased relative to actual planting practices and as presented in Table 7.

Conclusion
Tracking of daily crop growth and evapotranspiration within water districts is important for assessments of ET c within irrigated lands of the Lower Colorado River Basin.A satellite-based vegetation index model that accomplishes these tasks was developed, tested, and validated using remotely sensed Sentinel-2 NDVI data, Reclamation crop maps, local weather, and Google Earth Engine resources.The approach, 'Vegetation Index Model', constructs FAO-56 style crop growth stages from NDVI, and then estimates ET c .The aim is to improve and extend the methodology provided by the Evapotranspiration Report.Four crops, broccoli, cotton, wheat, and alfalfa-selected to represent short, middle and perennial seasons-were evaluated two ways: first by comparing models using Reclamation-verified field locations; and second, by validating models using ET c_act obtained from 13 eddy covariance sites.In the first way, ET c results from the Evapotranspiration report overestimated alfalfa and wheat by 21-25%, overestimated cotton by 6%, and underestimated broccoli ET c by 21%.The Vegetation Index Model showed ET c variability ranged 6-18% of median total ET c .In the second way, ET c results showed discrepancies ranging 1-14%, indicating that Vegetation Index Model results were valid.Results also indicated that the Vegetation Index Model is more representative of field-scale ET c than values provided by either the Evapotranspiration Report or by SIMS.By using remote sensing to adaptively estimate ET c for every field, ET c estimation bias can be reduced, and field-to-field variability quantified.
Results from the study show that revision of crop growth stages can be practical with available satellite data and software tools.Extension of the Vegetation Index Model to other irrigation districts can be direct if accurate weather data and crop maps are available.The Model would be more difficult to apply for locales with predominantly cloudy skies.The crop growth mapping approach can be aggregated to district-wide averages, and it is well suited for field-by-field applications.This latter option would be especially helpful for meeting Reclamation's obligations to provide accurate annual assessments.
Future work is needed to evaluate growth stages and ET c for crops not evaluated in this study, including lettuce, citrus, dates, and Bermuda grass.Also needed will be comparisons between the Vegetation Index Model and surface energy balance models, as provided by the OpenET Project (Volk et al., 2021;Volk et al., 2023).These comparisons will help resolve which ET mapping approaches-vegetation index, surface energy balance, or some combination of them-are most practical, accurate, and cost-effective.One avenue for further model improvement could lie with the use of active radar remote sensing data, as provided by Sentinel-1 and by the upcoming NISAR mission (nisar.jpl.nasa.gov).These sensors can discriminate cropped from non-cropped areas, detect early season irrigation events (e.g., Bazzi et al., 2022) and are unaffected by cloud cover, all attributes that would improve modeling of ET c immediately following planting and help extend use of the Vegetation Index Model into more humid climates.

Fig. 1 .
Fig. 1.Time series comparison of crop coefficients (Kc x 10, green), reference evapotranspiration (ET os , red), and crop evapotranspiration (ET c , black), illustrated for the Evapotranspiration Report small grains (wheat) model.Transitions between growth stages are delineated by vertical dotted lines.

Fig. 2 .
Fig. 2. Study area in Arizona (AZ) and California (CA), extending 100 km east-west and 40 km north-south.Mexico lies to the south and west.Reclamation-supplied field validated polygons (70,000 ha) are outlined in red, other fields are outlined in gray.Extent of the Yuma-Wellton evaluation area is outlined in blue, AZMET weather stations indicated by black symbols, and eddy covariance sites by yellow.

Fig. 5 .
Fig. 5. Partitioning alfalfa cuttings at Yuma.Observed NDVI values indicated with solid black symbols, NDVI minima by open circles, and NDVI maxima by triangles.The intersections of a seasonal trendline (red) with the NDVI trace facilitate discrimination of local minima and maxima and cutting events (vertical dashed lines).

Fig. 6 .
Fig.6.Sensitivity tests of crop growth stages (top row) and single crop coefficients, K c (bottom row) for broccoli, cotton, and wheat in 2019.
A.N.French et al.

Fig. 7 .
Fig. 7. NDVI traces for every ground-validated alfalfa, broccoli, cotton, and wheat site grown at Yuma in 2020 (gray).Example traces for each crop are shown in solid blue.Nominal Evapotranspiration Report K c curves are indicated by green lines.Year transition times denoted by dashed vertical lines.

Fig. 9 .
Fig. 9. Number of annual alfalfa cuttings at Yuma for 2018-2020 as resolved by the Vegetation Index Model.
evaluated 7 irrigation levels for one year of broccoli.For cotton, Anapalli et al. (2020) evaluated one 250 ha field over two years, Hunsaker and Bronson (2021) evaluated 8 fields over three years, and Ko et al. (2009) evaluated 6 fields over two years.Wheat K c values were studied on six fields for 2005-2008 by

Fig. 12 .
Fig. 12.Comparison of observed ET c_act (EC) vs. modeled daily ET c for example broccoli, wheat, cotton, and alfalfa sites.Results are shown for the Evapotranspiration Report (ER, green), SIMS (blue), and the Vegetation Index Model (VI, red).

Fig. 13 .
Fig.13.Modeled crop coefficient curve comparisons over Reclamation fields for broccoli, wheat, cotton, and alfalfa using the Evapotranspiration Report (K c ER, green), and Vegetation Index Model (K c VI, gray).Alfalfa curves are omitted for clarity.K c values corresponding to the eddy covariance validation sites(Table 8,  2018(Table 8,   -2020) )  are shown in red.For clarity, only one alfalfa site (WMIDD21-22b) is shown.Vertical line denotes January 1.
. Also reported are summaries by state and Mexico for consumptive use, bypassed water, banked and stored water, intentional surplus, and drought contingency contributions.The chief use for the Water Report is to document Decree entitlement compliance by users and, more recently, to document compliance with reduced allotments.The Water Report is scrutinized by users and entities in all seven Basin States.As the water supply is less than demand, every user, and especially the junior priority users, Metropolitan Water District of Southern California and the Central Arizona Water Conservation District, are particularly vigilant in reviewing all values in the Water Report.A companion report, initiated in the mid-1990's, entitled the 'Lower Colorado River Annual Summary of Evapotranspiration and Evaporation' (Evapotranspiration Report; available from https://www.usbr.gov/lc/region/g4000/4200Rpts/DecreeRpt/2015/2015.pdf and usbr.gov/lc/region/g4000/wtracct.html),adds details not provided in the Water Report.These include spatial analysis of agricultural and riparian vegetation water use, open water evaporation, crop type distributions, and estimated crop water use.The Evapotranspiration Report associates consumed water with users, entities, and districts.Most relevant to this study, it combines remote sensing observations of land cover type, weather data, and standardized crop evapotranspiration models, to estimate consumptive use.

Table 2
Vegetation Index Model implementation components.

Table 4
Evapotranspiration Report crop coefficients and growth stages for broccoli, cotton, wheat.
A.N.French et al.

Table 5
Evapotranspiration Report growth stages for alfalfa.Day of Year (DOY) are indicated for the first day of the cutting number.The first and last cuttings of a year overlap with the adjacent year so that 10 cuttings occur over 412 days.L END is set to zero days.K c INI and K c MID are 0.368 and 1.104, respectively.
Fig. 3. Comparing crop coefficients (green, K c values x10), reference evapotranspiration (red, ET os , mm/day) derived from Yuma Valley AZMET, and Sentinel-2 NDVI (blue, x10) patterns for four researched crops: wheat, cotton, broccoli, and alfalfa.Day of year annotations for broccoli and wheat are extended beyond nominal end of year (365) for plotting continuity.

Table 6
Quantile procedure for charting wheat, cotton, and broccoli stage length values.All steps are applied to each field individually.
Delbart et al., 2015, Cong, 2012)lant emergence and transition from partial to full-cover.The INI/DEV transition uses a 10% threshold (significantly earlier than the 20% threshold used by others;Delbart et al., 2015, Cong, 2012), marks the approximate start of crop development, and corresponds approximately to the 10% canopy cover threshold used in FAO-56(Allen et al., 1998, p 95).The beginning and end of the MID stage is selected by a 90% quantile, a value chosen to approximate the start of full canopy cover.With experience, optimal percentages for each crop might be identified.The transition days for INI/DEV, DEV/MID, and MID/END are found by linearly interpolating the 10% and 90% quantiles, which means that it is important to have locally smooth and non-replicated NDVI values.Accurate detection of the day of planting is a difficulty facing all NDVI approaches.It might correspond to an early season NDVI minimum, as used in

Table 7
Quantile procedure for charting alfalfa 'L' values.All steps are applied to each field individually.

Table 8
Eddy covariance site data used for ET c validation.

Table 9
Evapotranspiration Report and Vegetation Index Model estimated growth stages at Yuma, 2018-2020; where: Nnumber of analyzed sites, ER -Evapotranspiration Report, VI -Vegetation Index Model, Δ -difference in length of growth stages in days between VI and ER.

Table 10
Cutting durations at Yuma for Reclamation-verified sites 2018-2020.Variability is estimated by median absolute deviation (MAD).

Table 12
Day of planting estimation error at eddy covariance sites.Listed for each site are actual day of planting, its estimated value from NDVI minima, from FAO-56 INI values, and their corresponding errors.

Table 13
Evaluation of the efficacy of NDVI minima to select day of plantings at Yuma Reclamation-verified sites.Listed of the fraction of cases where the day of year at NDVI minima lie within the FAO-56 time selection window.

Table 15
Crop validation sites at Yuma, 2016-2022.Eddy covariance quality at each site is summarized by energy balance closure analysis, where regression results from comparing eddy covariance advective to radiative daily energy fluxes are reported by slope (b1) and R 2 .Listed are total season ET c _act from eddy covariance measurements (EC), followed by modeled ET c: Evapotranspiration Report (ER), Vegetation Index Model (VI), and SIMS.

Table 16
ET c model regression results for validation assessments of the Evapotranspiration Report (ER), the Vegetation Index Model (VI), and SIMS.