Integrating SEBAL with in-Field Crop Water Status Measurement for Precision Irrigation Applications—A Case Study

Abstract

The surface energy balance algorithm for land (SEBAL) has been demonstrated to provide accurate estimates of crop evapotranspiration (ET) and yield at different spatial scales even under highly heterogeneous conditions. However, validation of the SEBAL using in-field direct and indirect measurements of plant water status is a necessary step before deploying the algorithm as an irrigation scheduling tool. To this end, a study was conducted in a maize field located near the Venice Lagoon area in Italy. The experimental area was irrigated using a 274 m long variable rate irrigation (VRI) system with 25-m sections. Three irrigation management zones (IMZs; high, medium and low irrigation requirement zones) were defined combining soil texture and normalized difference vegetation index (NDVI) data. Soil moisture sensors were installed in the different IMZs and used to schedule irrigation. In addition, SEBAL-based actual evapotranspiration (ET_r) and biomass estimates were calculated throughout the season. VRI management allowed crop water demand to be matched, saving up to 42 mm (−16%) of water when compared to uniform irrigation rates. The high irrigation amounts applied during the growing season to avoid water stress resulted in no significant differences among the IMZs. SEBAL-based biomass estimates agreed with in-season measurements at 72, 105 and 112 days after planting (DAP; r² = 0.87). Seasonal ET matched the spatial variability observed in the measured yield map at harvest. Moreover, the SEBAL-derived yield map largely agreed with the measured yield map with relative errors of 0.3% among the IMZs and of 1% (0.21 t ha⁻¹) for the whole field. While the FAO method-based stress coefficient (K_s) never dropped below the optimum condition (K_s = 1) for all the IMZs and the uniform zone, SEBAL K_s was sensitive to changes in water status and remained below 1 during most of the growing season. Using SEBAL to capture the daily spatial variation in crop water needs and growth would enable the definition of transient, dynamic IMZs. This allows farmers to apply proper irrigation amounts increasing water use efficiency.

1. Introduction

Agriculture is a major and inefficient consumer of fresh water and the competition with other sectors for water is likely to increase exponentially in the next decades [1]. For this reason, a more reasonable and efficient use of water in agriculture is a crucial aspect in order to meet the food needs of the rapidly increasing world population [2,3].

Variable rate irrigation (VRI) has been demonstrated to provide better water use efficiency (WUE) while increasing crop yields [4,5,6]. VRI systems are today integrated with a wide variety of information in order to define proper rates for each irrigation management zone (IMZ), ranging from soil moisture sensing to thermal imagery [7,8,9,10].

Crop evapotranspiration (ET_c) is an important factor in estimating crop water demand and has been used to estimate yields and schedule irrigation [8,11]. Commonly, ET_c from agricultural areas is calculated by multiplying the reference evapotranspiration (ET₀; [12]) by a constant crop coefficient without considering the possible ET_c spatial variations in water or nutrient-limited zones ([13]).

Energy balance techniques [14,15,16,17] have been extensively used to map evapotranspiration at the local or regional scale through satellite-based image-processing tools such as the surface energy balance algorithms for land (SEBAL; [18]) and mapping evapotranspiration at high resolution with internalized calibration (METRIC™; [13]). These approaches have been tested under different soil and climate conditions in more than 30 countries worldwide [19,20,21,22] with typical field scale accuracy in ET estimates ranging between 85% for one single-day events up to 95% on a seasonal basis [20].

SEBAL estimates the evapotranspiration component converting the sensed satellite radiances (visible, near infrared and thermal-infrared) into surface parameters such as surface temperature, surface albedo, and normalized difference vegetation index (NDVI). The different components of the energy balance equation are then calculated using these surface parameters [13]. SEBAL’s main advantages are (i) limited input data required, (ii) it estimates ET without a priori knowledge on soil, crop and management conditions (no land use classification required), (iii) it is based on a physical concept, applicable for various climates and (iv) it can be applied at different spatial and temporal resolution [18].

While remote sensing-based energy balance techniques (e.g., SEBAL) have been extensively used to map evapotranspiration at different scales with interesting results [21,22], the effort in using remotely sensed data to evaluate water management practices and schedule irrigation has been limited. Some authors [23] argued the temporal and spatial resolution of currently available remote sensing data are not sufficient to use actual evapotranspiration estimates for on-farm research. Typical turnaround time for satellite image acquisition and processing is between one and three weeks [23], whereas ET maps should be delivered almost in real-time to correctly use the remotely sensed information for irrigation scheduling. Moreover, the spatial resolution of ET maps calculated applying remote sensing-based energy balance techniques frequently does not match with spatial resolution requested for VRI.

However, data fusion, sub-pixel extraction [23] and temporal integration techniques [24,25,26] could represent an alternative to improve both spatial and temporal resolution. The availability of new satellite platforms with higher temporal and spatial resolution is a need for a better deployment of these models. In addition, evapotranspiration estimates need to be validated by means of well-established methods, such as soil moisture sensing, shown to be a strong indicator of plant water status in order to define specific irrigation thresholds [10].

The objective of this study is to calculate SEBAL-based actual ET estimates in a pivot-irrigated field in North-Eastern Italy and compare the results to soil moisture data and crop stress coefficients (K_s). Accuracy of the SEBAL results is then assessed by comparing spatially distributed biomass and crop yield measurements.

2. Materials and Methods

2.1. Study Area and Plant Material

To validate SEBAL-based ET and biomass estimates, a study was conducted on a pivot-irrigated, 16-ha field located in Pettorazza Grimani, Italy (45°6′7.19″N, 12°0′45.51″E; Figure 1). Maize (Zea mais L.) seeds of PR32B10 (Pioneer Hi-Bred International, Inc., Johnston, IA, USA) variety were planted on April 13th, 2015 at 0.45 m inter-row spacing and seeding rate of 8 seeds m⁻². Prior to planting, cultivation operations included an autumn plowing at 30 cm depth and a seedbed preparation carried out by using a spring tine cultivator. Fertilization included a base-dressing of 65 kg N ha⁻¹ with manure, 92 kg P₂O₅ ha⁻¹, 111 kg K ha⁻¹ K₂O before the autumn plowing and a top-dressing of 184 kg N ha⁻¹ (urea). Maize was harvested on September 7th, 2015.

2.2. Soil Sampling and Management Zone Delineation

Irrigation management zones were defined integrating soil texture and NDVI data. A mixed-sampling scheme of the top soil layer (0–0.3 m) was followed based on a regular grid: 55 samples were collected at the nodes of 55 m × 65 m grids and 48 additional points at the nodes of 16 transects at 2.5, 15 and 30 m from each randomly chosen grid node, resulting in a total of 103 samples (Figure 2a; [27]). Soil texture was determined for all the 103 disturbed soil samples (Figure 2b,c) using a Malvern Mastersizer 2000 (Malvern Instruments Ltd, Great Malvern, UK) laser particle size analyzer (see Supplementary Materials for texture analysis method) and spatially interpolated using ordinary kriging. Maize spatial variability in NDVI was estimated at the peak of plant vigor (tasseling/silking growth stage) from previous seasons using LANDSAT 7 images of June 29th, 2009 and June 24th, 2013 [28]. LANDSAT 7 images were selected since they have the same spatial resolution as the VRI center pivot (30 m x 30 m pixels). Moreover, these images were selected since they represented the maize variability in NDVI during two extremely different growing seasons: A rainy (2009) and a dry (2013) year.

The investigated area was then divided into three IMZs based on the management zone analyst 1.0’s fuzziness performance index (FPI) and the normalized classification entropy (NCE) [29]. As discussed in [29], the optimal number of zones (clusters) was identified when each index had its minimum representing the least membership sharing (FPI) or greatest amount of organization (NCE). In addition, a fourth IMZ was defined and treated as the uniform irrigation zone.

2.3. Soil Moisture Assessment and Variable Rate Irrigation System Assessment

To monitor the soil moisture condition in the field, three 5TE soil moisture sensors (Decagon Devices Inc., Pullman, WA, USA) were installed in each IMZ at three different depths, 0.1, 0.3 and 0.5 m below the soil surface. Each array of soil moisture sensors was connected to an EM50G wireless cellular data logger (Decagon Devices Inc., Pullman, WA, USA). Before field installation, all sensors were calibrated in the laboratory (accuracy of ± 4%; [28]). The sensors collected soil moisture data every 30 min. Data were averaged for each array of sensors to get daily averages for the three different depths and then a unique average for each location. In the three IMZs, irrigation was scheduled when volumetric water content (VWC) approached 35% for zone 1, 32% for zone 2 and 28.5% for zone 3. The specific irrigation thresholds were defined based on soil hydraulic parameters for the different IMZs. Having the zones different soil characteristics, the thresholds were defined accordingly. In the Uniform IMZ, irrigation was applied at the maximum rate (38 mm) each time one of the other IMZs was irrigated.

Irrigation was supplied using a 270-m long, center pivot equipped with a Valley^® (Valley Irrigation, Valley, USA) variable rate irrigation system. Each sprinkler was controlled (on or off) by a pneumatically actuated flow control valve. The 96 sprinklers were grouped into 10 zones, each about 25 m in length. The uniformity of the irrigation system was assessed before the beginning of the growing season, following the American Society of Agricultural and Biological Engineers (ASABE) standards and the “catch can” layout proposed by [28,30,31]. The uniformity of application of the conventionally operated pivot (without VRI) was evaluated using the Hermann and Hein (CU_HH) uniformity coefficient [32] and the lower quarter distribution uniformity (DU_lq) coefficient [30]. The uniformity test resulted in high CU_HH and DU_ld (93.9 and 0.91, respectively). Moreover, measured collector volumes (observed) were compared to the predicted volumes (expected) [33] in order to test the performance of the VRI system. The application rates were set at 50%, 60%, 70%, 80% and 90% of the maximum pivot rate (30 mm). The catch cans were placed following two patterns, transect and arc-wise. Root mean square error (RMSE) of 1.7 and r² = 0.86 highlighted a good fit between expected and observed rates [28].

2.4. Evapotranspiration, Biomass and Yield Estimates Using the Surface Energy Balance Algorithm for Land (SEBAL)

The surface energy balance algorithm for land (SEBAL) determines the actual evapotranspiration (ET_r) of each pixel by solving the energy balance equation (Equation (1); [18]):

$$\mathsf{\lambda }\mathrm{E}={\mathrm{R}}_{\mathrm{n}}-{\mathrm{G}}_0-\mathrm{H},$$

(1)

where R_n is the net radiation, λE is the latent heat flux, G₀ is the soil heat flux and H is the sensible heat flux.

R_n is calculated from surface reflectance and surface temperature, the soil heat flux (G₀) as a function of R_n, surface temperature and NDVI. H is considered proportional to the ratio between the surface-air temperature difference (dT) and bulk aerodynamic resistance (r_ah) [18,24]. SEBAL computes the sensible heat flux H through a “calibration” procedure where dry (R_n = H + G₀, thus λE_dry ~ 0) and wet (R_n = G₀ + λE, H_wet ~ 0) pixels are selected and used to calculate the vertical temperature difference (ΔT) required to match the variation in H between the dry and wet pixels [26]. The wet pixel was identified as a very well watered crop pixel (low temperature, low albedo and high NDVI) and the dry pixel was searched between poorly vegetated or base soil areas with high temperature, high albedo and low NDVI.

Therefore, SEBAL calculates the evaporative fraction (Λ) considering the H and λE components ([18]; Equation (2)).

$$\mathsf{\Lambda }=\frac{\mathsf{\lambda }\mathrm{E}}{\left({\mathrm{R}}_{\mathrm{n}}-{\mathrm{G}}_0\right)}=\frac{{\mathrm{R}}_{\mathrm{n}}-{\mathrm{G}}_{\mathrm{o}}-\mathrm{H}}{{\mathrm{R}}_{\mathrm{n}}-{\mathrm{G}}_{\mathrm{o}}}$$

(2)

Moreover, the evaporative fraction (Λ) is assumed to be constant throughout the day. In fact, the difference between the Λ at the moment of image collection and the 24-h Λ average is considered as non-significant [24,26].

The latent heat flux λE is then calculated multiplying the Λ by the diurnal net radiation (R_n24h) as in Equation (3):

$$\mathsf{\lambda }\mathrm{E}={\mathsf{\Lambda }\times \mathrm{R}}_{\mathrm{n}}{}_{24h}.$$

(3)

Landsat 7 and 8 satellite images from April 19th, May 14th and 29th, June 6th and 30th, July 8th, 16th and 24th, August 9th and 17th, and September 2nd, 2015 were downloaded and automatically processed using a web-based platform called Tethys (Centrale Valutativa s.r.l.) to estimate the energy balance components. All dates had favorable clear-sky weather conditions. Visible bands (bands 1, 2, 3, 4, 5 and 7) were used for albedo (α), and vegetation index calculations. Thermal band (band 6) was used for surface temperature (T_s) and sensible heat (H). The spatial resolution was 30 m × 30 m on the visible bands and 60 m × 60 m on the thermal band. An atmospheric correction was performed to convert the top-of-atmosphere (TOA) radiance into surface reflectance using the second simulation of satellite signal in the solar spectrum algorithm (6S; [34]).

To overcome the lack of satellite images on a daily basis, actual evapotranspiration (ET_r) was calculated for each of the aforementioned satellite image acquisition dates. Then, ET_r on the date j was divided by the reference evapotranspiration (ET₀) on the same date j to calculate the fraction ETrFj. ET₀ was estimated using the FAO-56 method (Penman–Monteith equation) [35]. The fraction ETrFj = ET_r/ET₀ was then considered constant for the time period between j and the successive acquisition date j + 1 [24,26]. Seasonal actual evapotranspiration (ET_seasonal) was then calculated using Equation (4).

$$E{T}_{seasonal}={{\displaystyle \sum }}_j(ETr{F}_j\times {{\displaystyle \sum }}_{i=t_j}^{t_j+1}E{T}_0{}_i),$$

(4)

where t_j and t_j + 1 represents the time period between two consecutive image acquisition dates.

Maize aboveground biomass production (in kg ha⁻¹ day⁻¹) for each pixel of the satellite images was calculated from Tethys using the fraction of absorbed photosynthetically active radiation (fPAR) = absorbed PAR (APAR)/total PAR (PAR), the evaporative fraction (Λ) and light use efficiency (e’). fPAR was estimated from NDVI values as shown in Equation (5) [17]:

$$fPAR=-0.161+1.257\times NDVI.$$

(5)

The light use efficiency factor was set at 4 g/MJ [34]. Total biomass (BIO) production was converted into total above-ground biomass as total biomass (BIO) × 0.77 (Equation (6); [17]).

$$BIO=e^{\prime }\times \Lambda \times fPAR\times 0.84.$$

(6)

Maize yield was obtained multiplying the total aboveground biomass for each of the irrigation management zones by a specific harvest index (HI) factor equal to 0.5. Measured and SEBAL-derived yield averages were also calculated for each IMZ, uniform zone and for the rainfed portion of the field (top right corner). Moreover, differences between the measured and SEBAL-estimated yield were calculated for each pixel as in Equation (7):

$$Delta=\frac{SEBAL yield-observed yield}{observed yield}\times 100.$$

(7)

2.5. Crop Growth Data

To monitor crop growth during the season, biomass samples were collected on three sampling dates, at 72 (June 24th, 2015), 102 (July 24th, 2015) and 115 (August 6th, 2015) days after planting (DAP). On each sampling date, two 1 m² sections were destructively harvested for each IMZ. The biomass was then dried for 24 h in a forced-air oven to determine the total above-ground dry weight (DW) in kg ha⁻¹.

Yield data was collected at harvest using a New Holland CR8080 (New Holland Agriculture, Italy) combine harvester equipped with a yield monitor (grain mass flow and moisture sensor) and a differential global positioning system (DGPS). The harvester had a 7.5 m header and collected yield data every second. Raw data, expressed as 14% grain moisture content, was imported using SMS™ Advanced (AgLeader Technology, Ames, USA), cleaned to exclude field-edge effect and low yield values due to combine maneuvers and expressed as dry matter per hectare (DM ha⁻¹). The combine-derived data points were interpolated and resampled at 30 m × 30 m to match the spatial resolution of the SEBAL yield map.

2.6. SEBAL-Based and Soil Moisture Data-Based Water Stress Coefficient

Water stress coefficients (K_s) were calculated for each soil moisture sensor location using the moisture data collected at the three different depths. Water content at −33 KPa (field capacity) and −1500 KPa (permanent wilting point) were estimated from the soil characteristics of each IMZ using the “Rosetta” model [36]. Soil moisture-based K_s was then calculated on a daily basis using the FAO-56 method as follows [12]:

$$K_s=\frac{TAW-Dr}{RAW},$$

(8)

where TAW (mm) is the total available soil water in the root zone, Dr (mm) is root zone depletion and RAW is the fraction of the TAW that a crop can extract from the root zone before reaching water stress (RAW assumed 0.5 × TAW).

The SEBAL-based K_s (K_{s SEBAL}) was calculated by extracting the SEBAL actual K_c (K_{c SEBAL}) for each sensor location on a daily basis. Therefore, K_{s SEBAL} was divided by the basal crop coefficient to estimate the specific K_s (Equation (9)).

$$K_{s SEBAL}=\frac{K_c{}_{SEBAL}}{K_{cb}},$$

(9)

where K_cb is the basal crop coefficient estimated according to the FAO method [35] and K_csebal is the actual crop coefficient calculated from SEBAL.

K_{s SEBAL} is assumed to be equal to 1 when K_csebal ≥ K_cb. The K_s values were compared over the irrigation season to investigate a possible relationship between the two methods and the sensitivity of K_{c SEBAL} in capturing water stress-related changes in evapotranspiration.

2.7. Statistical Analyses

A regression analysis was used to evaluate SEBAL performances in estimating the biomass collected on the aforementioned sampling dates. The coefficient of determination was also calculated to describe the proportion of the total variance described by the model.

3. Results

3.1. Irrigation Management Zones (IMZs) Delineation

The field presented high variability in soil texture moving along the north-south axis with coefficient of variation of 25% for sand and 23% for clay. The sandy areas were located in the southern part while the north part of the field had high clay content (Figure 2). NDVI maps were strongly dependent on the season with the sandy areas showing high NDVI values in rainy seasons as in the 2009 map (Supplementary Materials Figure S1) and low values in dry seasons (2013 map, Supplementary Materials Figure S1) [28]). Soil texture (sand and clay; Figure 2) and NDVI data were integrated in order to delineate the IMZs using the aforementioned Management Zone Analyst software. The optimal number of management zones was selected based on the fuzziness performance index (FPI) and the normalized classification entropy (NCE) indices as shown in Figure 3 [29]. The cluster analysis identified a total of four zones but only the three IMZs under the pivot-irrigated area were considered for the experiment (Figure 4a): Zone 1, located in the north-eastern part of the field, was characterized by silty loam texture; Zone 2, located in the central part of the field, by loam texture, and Zone 3, in the south-western part of the field, by sandy-loam texture. The irrigation application map (Figure 4b) was developed matching the aforementioned management zones (zone 1, 2 and 3), and a fourth, uniform irrigated area to the resolution of the center pivot. The application map was then downloaded to the VRI controller and used for every irrigation event.

3.2. Soil Moisture Data and Irrigation Prescriptions

In 2015, rainfall totaled 272.8 mm between planting and harvest (Figure 5). The growing season was characterized by a dry April, July, August, September and a wet period between June and July. Total irrigation rates applied for each IMZ and the uniform zone are presented in

Table 1. Irrigation supplied during the 2015 growing season for all IMZs and the uniform zone.

Irrigation Event (DAP)	IMZ 1 (mm)	IMZ 2 (mm)	IMZ 3 (mm)	Uniform Zone (mm)
75	34	38	30	38
81	34	38	30	38
87	38	38	30	38
95	38	38	38	38
102	38	38	34	38
115	38	38	30	38
122	30	30	24	30
TOTAL	250	258	216	258

[28]. The high water demand during the growing reason resulted in a progressive decrease of the VWC below field capacity from mid-July (93 DAP) to the beginning of August, since the time required from the pivot for a complete rotation was the limiting factor (Figure 6).

As discussed before, irrigation was triggered when VWC approached 35% for zone 1, 32% for zone 2 and 28.5% for zone 3 (black dotted lines reported in Figure 6). Moreover, VRI allowed water saving of 8 mm (IMZ 1) and 42 mm (IMZ 3) throughout the season when compared to uniform irrigation rates, while IMZ 2 totaled the same amount of water supplied in the uniform zone (Table 1).

3.3. Measured and SEBAL-Based Data

Figure 7 shows the maize biomass values in t DM ha⁻¹ hand harvested for each IMZ on the three aforementioned sampling dates (blue bars). The high irrigation amounts applied during the growing season to avoid water stress in the experimental area resulted in no significant differences among the four irrigation zones. SEBAL-based biomass values extracted from the sensor locations at 72, 102 and 115 DAP resulted in a good agreement with observed data (r² = 0.87, Supplementary Materials Figure S2). However, as shown in Figure 7, SEBAL considerably underestimated the measured biomass at every sampling date. This underestimation in biomass could be related to the different spatial resolutions of the two systems: While biomass samples were collected on a 1 m² sections, SEBAL provided biomass estimates on 30 m pixels × 30 m pixels.

SEBAL-derived ET_r maps on June 2nd, July 2nd and August 6th (Figure 8) exhibited a marked temporal variation across the field (NDVI maps used to calculate ETr are shown in Supplementary Materials Figure S3). While early season ET_r maps as June 2nd, 2015 presented a limited variability in ET with values ranging between 3.4 and 4.3 mm/d (Figure 8a), late in the season the variability increased as in the case of August 6th, 2015 (Figure 8c) with ET values between 4.0 and 6.2 mm/d. In general, the central part of the field had greater ET_r values and the right-hand part had low ET_r.

This is confirmed from the seasonal (April 13th to September 7th) evapotranspiration map (ET_seasonal; Figure 9). Seasonal ET ranged between 471 mm to 554 mm with the lowest values in the top right part of the field and the highest values in the central area. ET_seasonal map largely matches the spatial variability of the measured yield map, with high seasonal ET associated with high yielding zones and low ET characterizing the low yielding zones (Figure 10a).

The SEBAL-based yield map (Figure 10b) had a similar spatial pattern when compared to the measured yield map. Measured yield values ranged between 5.2 to 13.2 t ha⁻¹, whereas SEBAL yield ranged between 7.5 to 11.1 t ha⁻¹. In both maps, low maize yields were generally observed in the right-hand part of the field, with the lowest yielding area at the top right corner. Moreover, a relative error map was generated as the difference between measured and SEBAL yield values on a pixel basis (Figure 10c). Almost all the pixels with relative errors higher than 25% are located outside the pivot-irrigated area in the rainfed portion of the field, indicating actual biomass or harvest index (HI) could have been overestimated. However, more than 75% of the map presented relative errors lower than 25% between measured and SEBAL yields with most of the pixels having errors below 15%. At the whole field, SEBAL underestimated the measured yield by less than 1 % (0.21 t DM ha⁻¹), which represents an excellent performance in agreement with [37].

In addition, measured and SEBAL yield were averaged at the IMZ scale, considering also the rainfed part of the field (top right corner;

Table 2. Irrigation supplied during the growing season for all the management zones and the uniform zone.

Irrigation Management Zone	SEBAL Above-Ground Biomass (t DM ha−1)	Harvest Index	Measured Yield (t DM ha−1)	SEBAL-Estimated Yield (t DM ha−1)	Difference (%)
IMZ 1	20.77	0.50	10.92	10.53	−3.6%
IMZ 2	21.23	0.50	10.47	10.55	0.8%
IMZ 3	20.66	0.50	10.53	10.33	−1.9%
Uniform	21.07	0.50	10.85	10.53	−2.9%
Non-irrigated Zone	20.02	0.50	9.39	9.99	6.4%

). Measured yield averages between 10.47 and 10.92 t ha⁻¹ in the irrigated areas and 9.39 t ha⁻¹ in the rainfed portion. Differences in yield were not significant between the VRI and the uniform irrigated areas, however, VRI resulted in water savings up to 42 mm as mentioned in 3.2. SEBAL was in a good agreement with measured data, slightly underestimating yield in IMZ 1, IMZ 2 and uniform zone, between −1.9% and −3.6%, and overestimating yield (0.8%) in IMZ 2. A 6.4% overestimation was observed in the rainfed zone. Overall, SEBAL underestimated the total yield by 0.3%, which represents again an excellent result.

3.4. SEBAL and Soil Moisture-Based K_s

The water stress coefficient (K_s) was used to estimate possible reduction in the potential evapotranspiration over the season. Daily K_s were calculated using the previously described FAO method and the SEBAL actual evapotranspiration. Figure 11 shows the FAO K_s (orange lines) and the SEBAL-based K_s (blue lines) values calculated for all the IMZs and the uniform area.

VWC never exceeded the FAO 50% TAW threshold in all the management zones (Figure 6) and in consequence FAO K_s was always equal to 1 [36]. Conversely, the SEBAL K_s showed changes over the growing season in all the management zones with values frequently lower than the optimum K_s = 1. Moreover, the trend in SEBAL K_s resulted similar in all the zones probably because of the high frequency of irrigation events allowed to maintain the crop well watered irrespectively of the soil variability.

4. Discussion

Weather conditions observed during the summer, did not allow us to fully exploit the VRI potentials. Indeed, exceptionally high temperatures and prolonged drought were observed starting before maize flowering to mid-August. These conditions smoothed the differences between IMZs in terms of water retention capabilities and, as a consequence, between irrigation volumes. Nevertheless, the VRI system coupled with the soil moisture sensors optimized the irrigation scheduling and maximized yields, allowing for fine tuning of the irrigation, particularly in IMZ 1 and IMZ 3. In turn, maize yield did not show any significant differences among the management zones, while water savings compared to the business-as-usual regime ranged from 8 mm (low irrigation requirement zone) to 40 mm (high irrigation requirement zone). In particular, it is worth pointing out that in the sandy area (i.e., IMZ 3), it was possible to obtain a reduction of irrigation volumes, still maintaining yield performances not different from the irrigated control treatment (uniform IMZ).

In the current study we supported the hypothesis that SEBAL estimated parameters such as actual evapotranspiration, crop coefficient and daily biomass accumulation could be used to understand the spatial and temporal variation in the maize water status. Moreover, whether the integration between SEBAL parameters and soil moisture information could represent a viable tool for scheduling irrigation and providing better crop yields while reducing the irrigation amounts.

SEBAL was able to describe 87% of in-field biomass variability, but average data were overestimated with respect to the actual values. These discrepancies could be related to the different spatial resolution between SEBAL and in situ operations. While the former calculated the biomass on a 30 m × 30 m grid, the hand harvested samples were collected on 1 m² plots. Improving the sampling pattern, collecting a more representative sample for each selected pixel, or using SEBAL parameters with higher spatial resolution could represent possible alternatives to increase the accuracy of estimates on both sides [17].

Conversely, SEBAL-estimated maize yield well correlated to the measured yield map in all the aforementioned management zones but also on a pixel basis (average error less than 1%; 0.21 t ha⁻¹). The relative error map (Figure 10c) reinforces the excellent performance of the algorithm as noted by [37]. Moreover, SEBAL exhibited the highest overestimation of maize yield in the rainfed area (+6.4%). The authors suggest this was probably related to uncertainties in harvest index (HI), which could be affected by plants stress conditions. In fact, the integration of a crop stress-related harvest index model could represent a way to provide better SEBAL-based yield estimates especially for zones under non optimal conditions, as rainfed areas, or when stress periods occur even in irrigated portions over the growing season [38].

ET_r and crop coefficient maps calculated by SEBAL, clearly highlight the potentialities of the model to optimize VRI. Unfortunately, it was not possible to compare the accuracy of SEBAL and FAO methods with respect to field ET measurements (e.g., lysimeters and flux towers). As discussed in many other studies, SEBAL resulted in underestimation of ET when compared to the FAO dual coefficient method [19,22,35]. In addition, FAO-56 tabulated K_c and K_cb has been shown to overestimate the observed crop coefficients and ET for humid regions up to 40% [39,40] due to the complexity of the interactions between several soil, weather and plant physiological factors. Moreover, weather conditions (e.g., cloud cover and net radiation) at the time of the satellite overpass and the selection of dry and wet pixels could strongly influence the final ET_c values calculated from SEBAL [41,42].

Anyway, SEBAL ET_r maps and K_s values showed a marked variability over the growing season even if the FAO method suggested no stress occurred (K_s = 1). While the widely-used FAO method does not account for spatial variability in crop evapotranspiration, the SEBAL ability to use pixel-specific estimates of plant water status represent a promising tool for the definition of transient irrigation management zones as the season moves, thus enabling dynamic variable rate irrigation.

5. Conclusions

In this study, the surface energy balance algorithm for land (SEBAL) was used to map crop evapotranspiration and biomass accumulation and estimate crop yield in a pivot-irrigated maize field. Moreover, SEBAL estimates of crop water stress were coupled with in-season soil moisture sensing data in order to define an irrigation scheduling strategy able to integrate these two approaches. Measured biomass and maize yield values showed a good agreement with SEBAL estimates, indicating the algorithm could be used to provide accurate estimates of crop growth and maize final production. The SEBAL sensitivity to the complex interactions between several soil, weather and plant physiological factors, allowed also to describe the spatiotemporal variability of crop water demands at the field scale. These results suggest that the integration of SEBAL with in-field soil moisture sensing could allow farmers to optimize irrigation scheduling for variable rate irrigation systems, overcoming the simplified FAO single or dual coefficients methods. However, major advancements are still needed for a successful use of SEBAL for irrigation purposes. The specific relationship between SEBAL K_s values and crop water status should be defined in order to identify the proper dynamic variable rate irrigation zones. SEBAL should be coupled with crop simulation models (e.g., DSSAT) to overcome its limited use during long cloudy periods. SEBAL applications should be based on satellites data with higher spatial and temporal resolutions.

Anyway, SEBAL adoption for precision agriculture irrigation will allow us to advance the irrigation science particularly in environments where water availability is scarce or the quality of irrigation water is poor.