Predicting Crop Evapotranspiration by Integrating Ground and Remote Sensors with Air Temperature Forecasts †

Water use efficiency in agriculture can be improved by implementing advisory systems that support on-farm irrigation scheduling, with reliable forecasts of the actual crop water requirements, where crop evapotranspiration (ETc) is the main component. The development of such advisory systems is highly dependent upon the availability of timely updated crop canopy parameters and weather forecasts several days in advance, at low operational costs. This study presents a methodology for forecasting ETc, based on crop parameters retrieved from multispectral images, data from ground weather sensors, and air temperature forecasts. Crop multispectral images are freely provided by recent satellite missions, with high spatial and temporal resolutions. Meteorological services broadcast air temperature forecasts with lead times of several days, at no subscription costs, and with high accuracy. The performance of the proposed methodology was applied at 18 sites of the Campania region in Italy, by exploiting the data of intensive field campaigns in the years 2014–2015. ETc measurements were forecast with a median bias of 0.2 mm, and a median root mean square error (RMSE) of 0.75 mm at the first day of forecast. At the 5th day of accumulated forecast, the median bias and RMSE become 1 mm and 2.75 mm, respectively. The forecast performances were proved to be as accurate and as precise as those provided with a complete set of forecasted weather variables.


Introduction
Optimal irrigation water management is one of the main challenges of precision agriculture, especially in open field crops, where farmers must deal with the uncertainty of both weather conditions and water availability, under climate change scenarios [1]. In many regions of the world, such as Mediterranean areas, where open field irrigation is practiced during dry seasons characterized by negligible rain contributions, irrigation water requirements are essentially driven by crop evapotranspiration.
Post-processing procedures based on a simple bias correction approach are particularly successful in improving the forecast performance of prognostic variables, such as temperature, at locations affected by systematic errors [15].
The second point (ii) descends from the consideration that real time forecasts of the weather variables required for ET c forecasting according to the Penman-Monteith equation, are available only by subscribing costly services with the weather forecast providers. Public (and free) numerical data of weather forecasts are generally restricted to air temperature [16], while other weather conditions are generally described in a qualitative manner (e.g., weather type), or according to predefined scale (e.g., Beaufort scale of wind force).
This study, as an extension of what had been anticipated in a previous conference paper [17], presents a detailed description and evaluation of a new methodology designed for forecasting the ET c of open field crops, by integrating crop parameters data retrieved from freely available satellite multispectral imagery, weather data provided by ground sensors, and temperature forecasts. The performance of the proposed methodology is evaluated in the Campania region, Southern Italy.

Study Area and Data from Ground Sensors
The study area is the Campania region located in Southern Italy, as shown in Figure 1. The region has a surface area of about 14,000 km 2 , and extends between the Tyrrhenian Sea and the Southern Apennines. Most of the region is characterized by dry-summer, subtropical climates, with hot and dry summers. The mean monthly temperature ranges from 20 • C to 30 • C in summer, and between 5 • C and 15 • C in winter. The maximum monthly precipitation values are recorded during November and December, the minimum values during July and August. The spatial variability of the weather conditions is highly affected by the complex orography of the region [18]. The 18 ground-based automatic weather stations (AWS) have been operating since 2008. These stations are managed by the Regional Meteorological Service, and obey the EUMETNET technical specifications [19]: each station is equipped with two redundant sensors to provide measurements with high accuracy and precision standards. The AWS sites were chosen to achieve a good representation of the weather variability within the region, including both the coastal areas and the  Irrigation in open fields does not start before April, frequently in June, and it lasts until the end of September. However, the actual time range of the irrigation season is determined by the weather fluctuations and specific agricultural practices: in the following, we assume an irrigation season for the analyzed study cases from June to September.
In this region, two sets of data have been combined in order to achieve a representative ensemble of cases for evaluating the performances of the ET c forecasted with the proposed methodology: • Remotely-assessed canopy parameters of maize crop fields belonging to two farms, called Improsta and Soffritti; • Weather observations and forecasts at 18 stations distributed across the Campania region.
The locations of the two farms and 18 weather stations are depicted in Figure 1.
The combined use of these two sets of data for evaluating the performance of the proposed methodology in forecasting ET c is formally legit, since the temporal pattern of the canopy parameters, even though these can change between the different locations due to changes in local weather conditions across the Campania region, are statistically independent of the local NWP errors, and thus of the simulated ET c forecasts.
The two farms Improsta and Soffritti are representative of two different topographic conditions: The Improsta farm is in a large floodplain; the Soffritti farm is in a hilly inland area ( Figure 1).
Two important field campaigns were conducted during the irrigation seasons from June to September 2014, in seven fields at the Improsta farm, and from June to September 2015 in six fields at the Soffritti farm. Ground-based measurements of the leaf area index (LAI), i.e., the one-sided green leaf area per unit ground surface area, were obtained with a LICOR LAI-2000 Plant Canopy Analyzer, and were used to validate the LAI retrieved from multispectral satellite imagery [6]. Ground validation of satellite LAI estimates are fundamental in areas with complex terrain features, which highly affect the quality of the multispectral satellite imagery [4].
The 18 ground-based automatic weather stations (AWS) have been operating since 2008. These stations are managed by the Regional Meteorological Service, and obey the EUMETNET technical specifications [19]: each station is equipped with two redundant sensors to provide measurements with high accuracy and precision standards. The AWS sites were chosen to achieve a good representation of the weather variability within the region, including both the coastal areas and the inland side of the region. The available data from these ground AWSs are air temperature and humidity at 2 m; global incoming solar radiation; wind speed at 10 m; barometric pressure and precipitation (that, for the scope of this study, has been neglected). All data have a time resolution of 10 minutes, and are then gathered in average daily values. These sites, characterized by significantly different terrain conditions (from coastal to internal areas and from floodplains to hilly and mountain areas), are also relevant for having a representative evaluation of the NWP performances applied to ET c forecasts. In fact, weather forecasts performances may vary even within such a small region as Campania, because NWP model resolutions cannot resolve the small-scale spatial variability of the weather conditions, especially in areas characterized by complex terrain features [13]. Table 1 presents the geographical coordinates of the AWS locations along with their elevations, ranging from 1 m a.s.l. to 848 m a.s.l., and the mean values of the weather data of interest for estimating ET c with Penman-Monteith equation, registered from June to September in the years 2010-2015. It is interesting to observe that wind speed exhibits a high variability across the region, with values ranging from 1.10 ms −1 to 4.42 ms −1 . This large variability of the wind speed reflects the impact of the local morphology. The highest wind speed values are recorded by those AWS (n. 5, 10 and 12) located close to the ridges of the Southern Apennines, exposed to winds blowing between the west and the eastern slopes. Average temperature ranges from 20.2 • C and 26 • C, with a variability mainly influenced by the distance from the Tyrrhenian seacoast and the elevation. The average relative humidity and solar radiation do not exhibit specific spatial patterns, being rather influenced by site-specific conditions.

Numerical Weather Forecasts
Forecasts of atmospheric pressure, solar radiation, wind speed at 10 m, relative humidity and temperature at 2 m were acquired from COSMO-LEPS outputs, with lead times from 1 to 5 days.
COSMO-LEPS is a limited area ensemble prediction system, developed within the Consortium for small-scale modeling (COSMO), and operated by the HydroMeteoClimate Regional Service of Emilia-Romagna, Italy (ARPA-SIMC) with 16 (that now are 20, since 2016) members, and a spatial resolution of 7.5 km.
The study in [10] firstly showed that COSMO-LEPS produces skillful and reliable forecasts up to 5 days of the weather variables relevant for ET c estimates in the Campania region. Moreover, it showed that the forecast performances do not significantly reduce by increasing the forecast lead time.
In this study, the numerical weather forecasts at the 18 AWS sites have been produced by a triangle-based, bi-linear interpolation method, which consists of interpolating the three numerical grid points closest to the examined AWS site. Then, the median value of the ensemble forecasts was computed for each AWS site. Since the COSMO-LEPS numerical outputs have a temporal resolution of three hours, the daily values of the forecasted variables have been calculated as the average values of the eight available forecasts produced for each day.
COSMO-LEPS outputs, as all numerical weather forecasts, exhibit systematic and non-systematic errors when compared with local weather observations. These forecast errors, originating from the inability of the NWP model to represent the atmospheric processes at the sub-grid scales, can be reduced by using different statistical post-processing techniques when the weather data observed at ground AWSs are available [20].
Here, a simple bias correction of weather forecasts has been performed to improve local forecasts. The bias correction technique consists in correcting the forecasts by the mean bias computed in the reference period [21]. The mean bias to be removed from the forecasted weather variables in the irrigation seasons 2014 and 2015 was estimated by comparing COSMO-LEPS forecasts and ground weather data at each AWSs during the irrigation seasons (June-September) of the years 2010-2013. Then, the bias-corrected temperature forecasts were obtained by subtracting these mean biases from the original temperature forecasts. Only temperature forecasts were bias corrected, since temperature was the only forecasted weather variable affected by systematic errors that could be effectively removed by the bias correction at each AWS site. The application of the bias correction to the other forecasted weather variables produced contrasting results across the study region. Although an average slight increase of the performance was observed across the region, at some locations the forecast performance decreased, due to the presence of non-systematic errors in areas where the topographic conditions limit the NWP model performances.

Remote Sensors and Multispectral Images for Retrieving Crop Parameters
Maize canopies at both the Improsta and Soffritti fields have been monitored by processing multispectral images, provided by remote sensors operating on DEIMOS-1 and LANDSAT 8 satellites. Images from Sentinel-2 satellites were not available for the period to which this study refers.
The crop canopy parameters used for estimating crop evapotranspiration with the Penman-Monteith equations are the albedo and the Leaf Area Index (LAI).
Albedo was used for estimating the fraction of incoming short-wave radiation leading to net radiation (R n ) in the Penman-Monteith equation, and it represents an estimate of the hemispheric and spectrally integrated surface albedo.
LAI was exploited for assessing the aerodynamic and surface resistances that appear in the Penman-Monteith equation [2].
LAI and albedo maps at Improsta Farm were retrieved by processing DEIMOS-1 multispectral satellite images with a spatial resolution of 22 m, and LANDSAT 8 multispectral images with a spatial resolution of 30 m, to obtain data with an average time frequency of 1 week, from June to September 2014. LAI and albedo maps at Soffritti Farm were derived from DEIMOS-1 and LANDSAT 8 multispectral imagery with a time frequency of 15 days, from June to September 2015.
DEIMOS-1 is a commercial tasking EO satellite, and is part of the Disaster Monitoring Constellation (DMC) (http://www.dmcii.com). The sensor records radiance in three spectral bands, corresponding to green (520-600 nm), red (630-690 nm) and near infrared (770-890 nm) parts of the electromagnetic spectrum at a ground sampling distance (GSD) of 22 m. DEIMOS-1 images were processed to a high degree of accuracy using an industry standard atmospheric correction algorithm (ATCOR-2) [22,23]. LANDSAT 8 OLI data at a GSD of 30 m were generated routinely by the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS), and obtained through the United States Geological Survey (USGS) Land Science Research and Development (LSRD) website; LEDAPS uses Dense Dark Vegetation (DDV) targets to estimate the aerosol optical thickness (AOT). The estimated AOT was used afterwards as input to the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer model.
Taking into account the minimal spectral resolution of the EO data normally available, the albedo (α) was calculated as a weighted total of the spectral surface reflectance (ρ bi ) for a given band bias resulting from atmospheric correction with broadband coefficients (ω bi ), reflecting the corresponding fraction of the solar irradiance in each sensor band [24]: LAI was retrieved from the multispectral by using the CLAIR model. The CLAIR model assumes an inverse exponential function between LAI and the weighted differential vegetation index (WDVI) [25,26]. The WDVI is calculated by means of the reflectance observed in the remote images in red and near-infrared bands.
LAI was then calculated as follows: where WDVI ∞ is the asymptotically limiting value for the WDVI, while b is an extinction and scattering coefficient. The two CLAIR model parameters (WDVI ∞ and the α* coefficient) can be identified by using chronologically paired LAI field measurements and reflectance values from satellite images [27]. In this study, LAI was computed by applying α* = 0.35, a value that was verified for this area from previous field campaigns [26].

FAO 56 Reference Evapotranspiration
The reference evapotranspiration is the evapotranspiration under standard conditions of a hypothetical grass reference crop, according to FAO Irrigation and drainage paper n. 56 [2]. It only depends on the weather data, and it was introduced to describe the evaporative demand of the atmosphere [2]. In this study, daily reference evapotranspiration was computed by applying daily weather data retrieved from ground sensors at the AWS sites. Hereinafter, we use the symbol ET 0-PM,obs to specify that reference evapotranspiration is computed with ground-based weather data, rather than numerical weather forecasts. ET 0-PM,obs (mm day −1 ) is computed as follows: where λ is the latent heat of vaporization, e s is the saturation vapor pressure (kPa), e a is the actual vapor pressure (kPa), ∆ is the slope of the vapor pressure curve (kPa • C −1 ), and γ is the psychometric constant (kPa • C −1 ). T is the daily mean air temperature at 2 m height ( • C). WS is the wind speed at 2 m height above ground (m s −1 ), that can be computed from the wind speed at 10 m above the ground by employing the logarithmic equation of the wind speed profile suggested by Allen et al. [2], that gives a multiplying factor of about 0.75. R n is the net radiation at the crop surface (MJ m −2 day −1 ). and G is the soil heat flux density (MJ m −2 day −1 ). The net radiation (R n ) was calculated as the difference between the incoming net shortwave radiation and the outgoing net long-wave radiation. As suggested by Allen et al. [2], for the reference crop, the incoming net shortwave radiation is a fraction of the incoming shortwave solar radiation (RS), by setting the albedo equal to 0.23. Daily mean air temperature (T) was computed as the average of daily maximum and minimum air temperature computed at the highest available temporal resolution of data. The actual vapor pressure was computed as a function of the mean air relative humidity (RH).

Hargreaves-Samani Equation for Estimating Reference Evapotranspiration
The Hargreaves-Samani equation was recommended as simplified method for estimating the reference evapotranspiration [2] for those sites where the availability of reliable weather data is questionable.
Since Hargreaves-Samani is an empirical equation, its prediction performance varies according to the climatic characteristics of the region where it is used. Thus, its local calibration is recommended to better estimate ET 0 [2]. Past studies suggested different techniques for calibrating the Hargreaves-Samani equation. The most common approach consists in calibrating a scaling coefficient [2,28].
Hargreaves-Samani referenced evapotranspiration is thus defined according to the following expression: In Equation (4), ET 0-HS is the Hargreaves-Samani reference evapotranspiration in (mm day −1 ), R a is the extraterrestrial radiation (MJ m −2 day −1 ), T max and T min are respectively the daily maximum and minimum temperature ( • C). The parameter k HS is the scaling factor to be calibrated.
Hereinafter, we use the symbol ET 0-HS,obs when Equation (4) is applied with air temperature retrieved from ground sensors, while ET 0-HS is computed with forecasted air temperature.
The scaling factor k HS is calibrated by minimizing the squared differences (ET 0-HS,obs − ET 0-PM,obs ) 2 summed over the entire calibration period. In this study, k HS was locally calibrated at each AWS by

Application of the Penman-Monteith Equation with Remotely Assessed Crop Parameters
Following the remote sensing-based Penman-Monteith direct approach (also known as the one-step approach, [3]), the Penman-Monteith equation was employed for estimating the crop potential evapotranspiration under standard conditions as follows: where ρ is the mean air density at constant pressure (kg m −3 ), c p is the specific heat of the air (MJ kg −1 • C −1 ), r a and r s are respectively the aerodynamic and surface resistances (s m −1 ).
The weather data in Equation (5) can be either derived by ground sensors or by weather forecasts. Herein after we use the symbol ET c-PM,obs to indicate the crop evapotranspiration estimated according to the Equation (5) with observed weather variables, reserving the ET c-PM when Equation (5) is applied with forecasted weather variables.
The aerodynamic resistance is defined as follows: where h c is the crop height, that was set equal to 0.4 m, as common practice. It has been evaluated that the sensitivity of ET c estimates to crop height in the range 0.1-0.6 m is negligible [25].
The surface resistance was computed as follows: where LAI was estimated from remotely-sensed crop images. The fraction of incoming short-wave radiation contributing to the net radiation R n was computed by multiplying the incoming short-wave radiation by (1-α), where α is albedo obtained by processing the satellite multispectral images.
The Albedo and LAI in Equations (5)- (7) were updated, as a new multispectral crop image has been processed, and they were kept constant until the next satellite observation becomes available, as done in most operational advisory services [6]. Crop growth models could be applied to predict the time variation of crop parameters between consecutive satellite passages, but this is beyond the scope of this paper.

Crop Coefficient
The crop coefficient, K c , is a dimensionless coefficient that is defined as the ratio between the evapotranspiration of the examined field crop under standard conditions and the reference crop evapotranspiration.
Here, we were interested in analytically computing the crop coefficient K c as follows: where ET 0-PM,obs and ET c-PM,obs , respectively, refer to the reference and the crop evapotranspiration computed with the Penman-Monteith equation with ground-measured weather data. ET c-PM,obs was computed with Equations (5)- (7), with the crop canopy parameters LAI and albedo estimated from multispectral satellite images. The crop coefficient, as defined by Equation (8), depends on both crop parameters and weather data. In the proposed methodology, crop parameters in K c are updated as a new crop canopy image becomes available (almost weekly), while weather data are daily updated based on ground weather sensors.

Forecasting Crop Potential Evapotranspiration with Hargreaves-Samani Equation
In the proposed methodology, ET c forecasts are estimated for lead times ∆t from 1 to 5 days by employing the Hargreaves-Samani equation (Equation (4)) with the forecasted bias-corrected temperature (instead of the observed temperature), multiplied by K c updated at the time t 0 of forecast delivery: ET c- For the sake of comparison, crop potential evapotranspiration has been also computed with the one-step approach by means of the Penman-Monteith ET c-PM applied with the forecast variables, as done by Chirico et al. [6]. Figure 2 shows a flowchart that provides an overview of the proposed methodology. The input data sources are on the top of the diagram: (i) observed weather data, (ii) multispectral images and (iii) weather forecasts Our methodology intends to provide ET c-HS forecast, by exploiting only temperature forecasts. The flowchart also illustrates the methodology for ET c forecasting based on Penman-Monteith, thus exploiting a complete set of weather variables to be acquired from meteorological services.

Forecasting Crop Potential Evapotranspiration with Hargreaves-Samani Equation
In the proposed methodology, ETc forecasts are estimated for lead times Δt from 1 to 5 days by employing the Hargreaves-Samani equation (Equation (4)) with the forecasted bias-corrected temperature (instead of the observed temperature), multiplied by Kc updated at the time t0 of forecast delivery: For the sake of comparison, crop potential evapotranspiration has been also computed with the one-step approach by means of the Penman-Monteith ETc-PM applied with the forecast variables, as done by Chirico et al. [6]. Figure 2 shows a flowchart that provides an overview of the proposed methodology. The input data sources are on the top of the diagram: (i) observed weather data, (ii) multispectral images and (iii) weather forecasts Our methodology intends to provide ETc-HS forecast, by exploiting only temperature forecasts. The flowchart also illustrates the methodology for ETc forecasting based on Penman-Monteith, thus exploiting a complete set of weather variables to be acquired from meteorological services.

Evaluation of Forecast Performances
Forecasts performances of the individual weather variables, as well as of ETc-PM and ETc-HS forecasts, have been evaluated by comparing them with the corresponding variables retrieved from

Evaluation of Forecast Performances
Forecasts performances of the individual weather variables, as well as of ET c-PM and ET c-HS forecasts, have been evaluated by comparing them with the corresponding variables retrieved from ground sensors. For what concerns the evapotranspiration, ET c-PM and ET c-HS forecasts are compared with ET c-PM,obs , hereinafter also denoted as "best estimate".
Statistical performance indices were computed for all lead times by comparing the median value of the ensemble of the forecasts on the generic i-th day, with the best estimate at the same day. The first index is the BIAS, which was used as an indicator of the accuracy of the forecasts: where X i f and X i obs , respectively, are the variables forecasted and observed at the i-th day, while n denotes the number of examined days, that in this study is 122 (from June, 1 to September, 30).
The second performance indicator is the root mean square error, RMSE, which gives insight into both the accuracy and precision of the forecasts. RMSE is calculated as follows: where the symbols are the same as those used in Equation (10). Forecast performances were evaluated at different lead times, i.e., from one day up to the five days after the weather forecasts were issued. The forecast performance of the evapotranspiration was evaluated by computing the corresponding values accumulated up to one, three and five days after the weather forecasts were issued.

Weather Forecast Performances
The RMSE of the forecasted weather variables was computed to better understand the performance of the evapotranspiration predicted with the Hargreaves-Samani and Penman-Monteith equations. Temperature forecasts were corrected with a simple bias correction technique to remove the systematic bias. The systematic bias was estimated in the months between June and September of the years 2010-2013 (reference period). The average (among all sites and lead times) mean bias of the raw temperature forecasts in the reference period is about −0.68 • C. The mean bias is always negative by varying site and lead time, thus denoting a systematic underestimation of the observed temperature by the forecast outputs of COSMO-LEPS. The underestimation was corrected by subtracting the mean bias from the raw temperature forecasts at each AWSs and for each lead time. This correction led to a significant reduction of the mean bias (almost equal to zero), related to the corrected temperature forecasts and of the RMSE from an average value of 1.70 • C to 1.16 • C (as better detailed for each site, lead time and year in Figure 3). The bias correction was not applied to the other variables, since it produced contrasting results across the region. The colored grids in Figure 3 depict the RMSE values calculated at the 18 AWS sites (along the rows) for five forecast lead times (from 1 to 5 days). The RMSE value is represented according to the color scale on the side: blue corresponds to the minimum RMSE; dark red to the maximum RMSE.
The forecast performance considerably changes across the region, but also from year 2014 to year 2015. Summer 2015 was characterized by anomalous unstable meteorological conditions, with colder temperature and less solar radiation. Temperature (T) RMSE is higher at sites like AWS n. 5-7 with irregular topography and at AWS n. 15 close to the coastline. Temperature RMSE is generally higher in computed with all variables produced by the weather forecasts (i.e., no observed weather variables were used for the RMSE assessment), in the year 2014 (light blue) and year 2015 (dark blue). The other couples of box plots indicate the RMSE of ETc-PM after substituting one forecast variable with the corresponding ground sensor measurement, while using the weather forecasts for the other variables.
The largest RMSE reduction is observed when the ground sensor data are used to substitute the relative humidity, i.e., the largest impact is due to errors in the forecasted relative humidity, which affects the vapor pressure deficit of the air relevant for the aerodynamic component of the Penman-Monteith equation.  Values are considered as outliers if they are larger than p75 + 1.5(p75 -p25) or smaller than p25 -1.5(p75 -p25).

Calibration and Validation of the Hargreaves-Samani Equation
The uncalibrated Hargreaves-Samani equation overestimates the reference evapotranspiration at most AWS sites, i.e., ET0-HS,obs computed with Equation (4) for kHS = 1 tends to be larger than ET0-PM,obs. As a result (see Table 2), the calibrated kHS is smaller than 1 at all AWS sites, except for AWS n. 5 and 10, where wind speed is generally very high, and thus the aerodynamic term (neglected by the Hargreaves-Samani equation) is more important. The average kHS is 0.85, ranging from 0.64 to 1.14. RMSE in the calibration period ranges from 0.51 to 1.04 mm/day, with an average value of 0.69   The largest RMSE reduction is observed when the ground sensor data are used to substitute the relative humidity, i.e., the largest impact is due to errors in the forecasted relative humidity, which affects the vapor pressure deficit of the air relevant for the aerodynamic component of the Penman-Monteith equation.

Calibration and Validation of the Hargreaves-Samani Equation
The uncalibrated Hargreaves-Samani equation overestimates the reference evapotranspiration at most AWS sites, i.e., ET0-HS,obs computed with Equation (4) for kHS = 1 tends to be larger than ET0-PM,obs. As a result (see Table 2), the calibrated kHS is smaller than 1 at all AWS sites, except for AWS n. 5 and 10, where wind speed is generally very high, and thus the aerodynamic term (neglected by the Hargreaves-Samani equation) is more important. The average kHS is 0.85, ranging from 0.64 to 1.14. The largest RMSE reduction is observed when the ground sensor data are used to substitute the relative humidity, i.e., the largest impact is due to errors in the forecasted relative humidity, which affects the vapor pressure deficit of the air relevant for the aerodynamic component of the Penman-Monteith equation.

Calibration and Validation of the Hargreaves-Samani Equation
The uncalibrated Hargreaves-Samani equation overestimates the reference evapotranspiration at most AWS sites, i.e., ET 0-HS,obs computed with Equation (4) for k HS = 1 tends to be larger than ET 0-PM,obs. As a result (see Table 2), the calibrated k HS is smaller than 1 at all AWS sites, except for AWS n. 5 and 10, where wind speed is generally very high, and thus the aerodynamic term (neglected by the Hargreaves-Samani equation) is more important. The average k HS is 0.85, ranging from 0.64 to 1.14. RMSE in the calibration period ranges from 0.51 to 1.04 mm/day, with an average value of 0.69 mm/day, while from 0.54 to 1.25 mm/day in the validation period, with an average value of 0.78 mm/h. Similar results were achieved by previous studies [29].  mm/day, while from 0.54 to 1.25 mm/day in the validation period, with an average value of 0.78 mm/h. Similar results were achieved by previous studies [29]. Figure 5 shows the scatter plot of the ET0-PM,obs versus ET0-HS,obs before and after the calibration applied to the validation dataset, of all AWSs in year 2014 and 2015. The scatter plots clearly show that the calibration of the kHS coefficient reduces both the bias and the spread of the predicted ET0-HS,obs compared with ET0-PM,obs.

Remotely Assessed Crop Coefficients
LAI was the crop canopy parameter that mostly influenced the temporal pattern of the crop coefficient computed according to Equation (8). Figure 6 shows the box plots of the LAI retrieved at the reference maize fields. The variability of the LAI observed at each date reflects the spatial variability of maize growth and biomass accumulation among the different selected fields. Maize at

Remotely Assessed Crop Coefficients
LAI was the crop canopy parameter that mostly influenced the temporal pattern of the crop coefficient computed according to Equation (8). Figure 6 shows the box plots of the LAI retrieved at the reference maize fields. The variability of the LAI observed at each date reflects the spatial variability of maize growth and biomass accumulation among the different selected fields. Maize at Improsta is representative of the expected temporal pattern of a maize canopy in favorable weather conditions in Southern Italy, reaching the maximum value (5) at flowering at the beginning of August. LAI patterns observed at Soffritti were affected by the anomalous weather of the summer in the year 2015, which determined less vigorous development of the maize crop. These two data sets of Improsta and Soffritti were selected for evaluating the suggested procedure, as they are representative of two extremely different crop conditions, and because they were favorably validated with field LAI observations [6].
Albedo also changes along with the crop development, from 0.15 in the early stage of development to 0.20 at the maximum development, but with a less pronounced variability among the different sites within the same farm and years. Figure 7 shows an example of the temporal pattern of the crop coefficient Kc at the two farms computed with weather data observed at AWS n. 16. The gray band reflects the effect of the LAI variability at each farm, as depicted by the box plots in Figure 6. The crop coefficient in the early stage of development is around 0.5. In the mid-season it reaches 1.2 at Improsta and 1.1 at Soffritti.

Forecasts of Crop Evapotranspiration
The performance of the ETc forecasts was evaluated with an ensemble data obtained by combining observed and forecasted weather data at 18 sites, and the ensemble of crop canopy collected from multispectral images in year 2014 and 2015. These two data sets of Improsta and Soffritti were selected for evaluating the suggested procedure, as they are representative of two extremely different crop conditions, and because they were favorably validated with field LAI observations [6].
Albedo also changes along with the crop development, from 0.15 in the early stage of development to 0.20 at the maximum development, but with a less pronounced variability among the different sites within the same farm and years. Figure 7 shows an example of the temporal pattern of the crop coefficient K c at the two farms computed with weather data observed at AWS n. 16. The gray band reflects the effect of the LAI variability at each farm, as depicted by the box plots in Figure 6. The crop coefficient in the early stage of development is around 0.5. In the mid-season it reaches 1.2 at Improsta and 1.1 at Soffritti. These two data sets of Improsta and Soffritti were selected for evaluating the suggested procedure, as they are representative of two extremely different crop conditions, and because they were favorably validated with field LAI observations [6].
Albedo also changes along with the crop development, from 0.15 in the early stage of development to 0.20 at the maximum development, but with a less pronounced variability among the different sites within the same farm and years. Figure 7 shows an example of the temporal pattern of the crop coefficient Kc at the two farms computed with weather data observed at AWS n. 16. The gray band reflects the effect of the LAI variability at each farm, as depicted by the box plots in Figure 6. The crop coefficient in the early stage of development is around 0.5. In the mid-season it reaches 1.2 at Improsta and 1.1 at Soffritti.

Forecasts of Crop Evapotranspiration
The performance of the ETc forecasts was evaluated with an ensemble data obtained by combining observed and forecasted weather data at 18 sites, and the ensemble of crop canopy collected from multispectral images in year 2014 and 2015. Two methods were compared:

Forecasts of Crop Evapotranspiration
The performance of the ET c forecasts was evaluated with an ensemble data obtained by combining observed and forecasted weather data at 18 sites, and the ensemble of crop canopy collected from multispectral images in year 2014 and 2015.
Two methods were compared: • ET c-PM , crop evapotranspiration forecasted by exploiting NWP forecasts for the entire set of weather variables required by the Penman-Monteith equation; • ET c-HS , crop evapotranspiration forecasted by exploiting NWP forecasts only for the air temperature, according to the proposed methodology.
The benchmark is represented by (ET c-PM,obs ), i.e., the crop evapotranspiration estimated with Equations (5)-(7) by using weather input data retrieved from ground sensors data according to the one-step procedure [3,4].
The forecast performances were evaluated by analyzing ET c accumulated in 1, 3 and 5 days after the forecasts were issued, since the interest was to evaluate the error in the estimated accumulated crop water requirement, which is relevant for scheduling the irrigation with different time intervals. The benchmark is represented by (ETc-PM,obs), i.e., the crop evapotranspiration estimated with Equations (5)-(7) by using weather input data retrieved from ground sensors data according to the one-step procedure [3,4]. The forecast performances were evaluated by analyzing ETc accumulated in 1, 3 and 5 days after the forecasts were issued, since the interest was to evaluate the error in the estimated accumulated crop water requirement, which is relevant for scheduling the irrigation with different time intervals.  Values are considered as outliers if they arelarger than p75 + 1.5(p75 -p25) or smaller than p25 -1.5(p75 -p25).
The accumulated forecasts underestimate ETc in 2014, while they overestimate ETc in 2015. In absolute terms, the median BIAS is around 0.2 for the first day of forecast, and accumulates up to 1 mm at the 5 th day of forecast (i.e., almost 0.2 mm/day). ETc-HS BIAS is slightly smaller than ETc-PM in absolute terms. More interestingly, the spread of the ETc-PM BIAS is larger than the ETc-HS BIAS, especially in year 2014, revealing the impact of the larger forecast uncertainty due to the larger The accumulated forecasts underestimate ET c in 2014, while they overestimate ET c in 2015. In absolute terms, the median BIAS is around 0.2 for the first day of forecast, and accumulates up to 1 mm at the 5 th day of forecast (i.e., almost 0.2 mm/day). ET c-HS BIAS is slightly smaller than ET c-PM in absolute terms. More interestingly, the spread of the ET c-PM BIAS is larger than the ET c-HS BIAS, especially in year 2014, revealing the impact of the larger forecast uncertainty due to the larger number of weather variables involved. The results in year 2015 are affected by the anomalous summer weather conditions, with strong atmospheric instabilities generating thunderstorms and strong winds in the afternoon [30].
Median RMSE is around 0.75 mm for the first day of forecast, and increases up to around 2.75 mm at the 5 th day of the accumulated forecast. The RMSE in 2015 is generally larger than in 2014, due to the strong weather instabilities in the summer of 2015, as mentioned above, which were more challenging to be predicted by the NWP model. In all cases, the ET c-PM median RMSE is larger than ET c-HS median RMSE. As for the BIAS, in the year 2014, ET c-PM exhibits an RMSE spread much larger than ET c-HS , as a result of the larger number of weather forecast variables involved in ET c-PM . The two spreads are comparable in the year 2015, since ET c-HS predictions were more exposed to the unstable weather conditions of the summer of 2015.
Overall, the performance of ET c-HS forecasts, based on temperature forecasts only, are comparable or even better than ET c-PM forecasts, which instead require a larger set of forecast weather variables.

Discussion and Conclusions
Irrigation advisory services should be designed to meet the requirements of three different user categories [5]: farmers, water managers responsible for irrigation water allocation at the district level, and water management authorities. The first two categories require irrigation advices based upon forecasts of crop water requirements predicted days ahead, in order to allow the optimization of the water irrigation allocation and irrigation schedules, accounting for the actual crop water needs. In many irrigated areas, such as in the Southern Mediterranean countries, forecasting crop water requirements essentially means forecasting crop evapotranspiration, since rainfall and ground water contributions are negligible. A forecast horizon of at least 5-7 days is desirable from an operational perspective [4].
To this scope, two conditions are required: • Crop canopy parameters need to be updated with relatively high frequency; • Reliable weather forecasts should be available, possibly locally bias-corrected by means of ground weather observations.
These two conditions are met by the procedure proposed in this study, with the advantage of using just temperature as the forecast variable, as the input of the Hargreaves-Samani equation, instead of the entire set of weather variables required by the Penman-Monteith equation.
The recent ESA Sentinel-2 mission [31] opened the possibility to estimate crop parameters every 5 days with a resolution of 10 m with free multispectral images. The temporal evolution of the crop parameters within this time window is generally negligible, and thus up-to-date crop parameters can be assumed to be invariant within forecast horizons of 5-7 days [6,32]. Since this study refers to those years 2014-2015, it was not possible to use Sentinel-2 imagery for crop parameter estimates. However, by combining two different satellite imageries provided by remote sensors operating on DEIMOS-1 and LANDSAT 8 satellites, a temporal frequency of 7-15 days for updating data could be achieved. Moreover, the spectral images of these two satellites provided an effective reconstruction of the temporal pattern of the crop parameters, as confirmed by the field campaigns conducted in two experimental sites in Campania.
In this research, we used COSMO-LEPS forecasts, a Limited Area Model operating in Italy, and providing weather data at high spatial and temporal resolutions, with forecast horizon of 5-7 days. This category of NWP models is characterized by high reliability, but their forecasts are distributed at high costs of subscription [10]. We showed that these costs could be reduced by using a simplified ET c estimation method, which employs only temperature forecasts.
There are a few reasons supporting the advisability of developing forecasting procedures that need only air temperature forecasts, rather than the entire set of weather variables affecting ET c processes. The forecast performances are not homogenous, being highly dependent on the local terrain features. Air temperature is better predicted than other meteorological variables, since it is less affected by sub-grid variabilities.
Post-processing techniques for correcting weather forecasts are more effective for air temperature than for other weather variables [20]. Differently from other weather variables, meteorological services freely broadcast numerical data of air temperature forecasts, by publishing them on dedicated web platforms.
The methodology presented in this paper has been designed for those sites where reliable ground weather stations are available and the Hargreaves-Samani equation works reasonably well.
The procedure exploits the remote sensing-based Penman-Monteith direct approach (also known as one-step approach, [3]) for assessing analytically the current ET c and crop coefficient. The crop coefficient approach is instead exploited for forecasting ET c in combination with the Hargreaves-Samani equation.
The results at 18 sites across the Campania region have proven that the suggested procedure produces results that are as accurate and precise as those provided with a complete set of forecasted weather variables, with better performances at sites and in weather conditions where the radiative component of the ET c equation is dominant over the aerodynamic term. In fact, in windy areas, the prediction performances of the Hargreaves-Samani equation were worse. In examined cases, for example, the highest root mean square errors were obtained at AWS sites n. 5-10-12, characterized by average wind speed larger than 4 m/s. The prediction performances in summer 2015 were clearly worse than in 2014 for a lead time greater than 1 day, since summer 2015 was characterized by anomalous thunderstorms. This circumstance mainly worsened the performance of temperature forecasts with lead time greater than 1 day, due to the inability of NWP models to predict the convective storms and to describe their small-scale variability. As a consequence, a simplified method based only on temperature forecasts for estimating ET c resulted in being less accurate in predicting the evaporative demand of the atmosphere at sites affected by those thunderstorms. However, even in 2015, the proposed procedure produced prediction errors comparable with those obtained by implementing the Penman-Monteith equation.
In this perspective, possible improvements could be achieved by adopting more advanced correction procedures of both Hargreaves-Samani and of the forecasted air temperature. In this study, a simple bias correction was adopted for correcting air temperature forecast errors. However, more sophisticated procedures could be implemented, accounting for the non-stationarity of the forecasting error [20,33]. Similarly, the Hargreaves-Samani equation was calibrated by applying a scaling factor. However, other procedures could be also tested, by involving other constants of the Hargreaves-Samani equation into the calibration process [34,35].