A Weekly Indicator of Surface Moisture Status from Satellite Data for Operational Monitoring of Crop Conditions

The triangle method has been applied to derive a weekly indicator of evaporative fraction on vegetated areas in a temperate region in Northern Italy. Daily MODIS Aqua Land Surface Temperature (MYD11A1) data has been combined with air temperature maps and 8-day composite MODIS NDVI (MOD13Q1/MYD13Q1) data to estimate the Evaporative Fraction (EF) at 1 km resolution, on a daily basis. Measurements at two eddy covariance towers located within the study area have been exploited to assess the reliability of satellite based EF estimations as well as the robustness of input data. Weekly syntheses of the daily EF indicator (EFw) were then derived at regional scale for the years 2010, 2011 and 2012 as a proxy of overall surface moisture condition. EFw showed a temporal behavior consistent with growing cycles and agro-practices of the main crops cultivated in the study area (rice, forages and corn). Comparison with official regional corn yield data showed that variations in EFw cumulated over summer are related with crop production shortages induced by water scarcity. These results suggest that weekly-averaged EF estimated from MODIS data is sensible to water stress conditions and can be used as an indicator of crops’ moisture conditions at agronomical district level. Advantages and disadvantages of the proposed approach to provide information useful to issue operational near real time bulletins on crop conditions at regional scale are discussed.


Introduction
Early warning systems able to exploit indicators of cropland status are essential sources of information to improve yield forecasts and resource management [1]. Sustainable management at regional scale of water resources relies on crop water needs information, especially over scenarios of climate change in anthropized regions where water is over exploited and shared among multiple usages.
Since more than 80% of the water applied to agricultural lands can be consumed by evapotranspiration (ET) [2], monitoring temporal and spatial ET variations is necessary. Traditional methods for direct or indirect in situ measurement of ET, such as lysimeters, sap flow meters and eddy covariance towers, provide accurate point and/or parcel-scale values but spatially wind speed, vapor pressure, and surface available energy). According to the authors, this allows reducing the underestimation of EF typically obtained using the "constant EF" assumption [50].
Finally, it is important to remind that, being based on optical RS data, EF time series often present gaps due to cloud cover. Therefore, if complete and continuous (e.g., daily) EF estimates are needed (e.g., for assimilation in hydrological models) appropriate gap filling approaches should be implemented. For example, Rasmussen et al. [27], developed a specific solution to deal with this problem, allowing to obtain a complete daily coverage of EF at regional scale.

Study Contribution and Objective
The previous review highlighted the maturity of surface moisture status assessment from satellite data, and that the triangle method can be applied operationally for a study area with specific characteristics, if performing a proper data analysis. In this framework, this paper aims to develop an operational approach, based on the triangular method, to derive a weekly spatially explicit indicator of crop surface moisture conditions at regional scale. The implemented solution produces weekly EF map over the agricultural areas of Lombardy, Northern Italy, covering more than 1,200,000 ha, starting from moderate resolution daily EO data from MODIS sensor. EF estimates were produced for the years 2010 to 2012 and locally validated with in-situ eddy-covariance data. The utility and robustness of the proposed indicator for operational monitoring of crop conditions at regional scale was assessed by comparing seasonal synthesis of EF and NDVI with official yield data.

Study Area Description
The study was focused on agricultural areas of Lombardy, in Northern Italy. The AOI is a part of Po valley and spans about 12,000 km 2 . Topography of the study area is mostly flat, with altitude slowly degrading from about 150 m in the North West to about 20 m in the Eastern area. It is intensively cultivated with summer crops (i.e., corn (Zea mays L.), forage and rice (Oryza sativa L.)) and winter wheat (Triticum L.), which cover about 68% and 8% of the agricultural land, respectively. It is divided in 45 agronomical districts, showing marked differences in the main cultivated crops (Figure 1). Finally, it is important to remind that, being based on optical RS data, EF time series often present gaps due to cloud cover. Therefore, if complete and continuous (e.g., daily) EF estimates are needed (e.g., for assimilation in hydrological models) appropriate gap filling approaches should be implemented. For example, Rasmussen et al. [27], developed a specific solution to deal with this problem, allowing to obtain a complete daily coverage of EF at regional scale.

Study Contribution and Objective
The previous review highlighted the maturity of surface moisture status assessment from satellite data, and that the triangle method can be applied operationally for a study area with specific characteristics, if performing a proper data analysis. In this framework, this paper aims to develop an operational approach, based on the triangular method, to derive a weekly spatially explicit indicator of crop surface moisture conditions at regional scale. The implemented solution produces weekly EF map over the agricultural areas of Lombardy, Northern Italy, covering more than 1,200,000 ha, starting from moderate resolution daily EO data from MODIS sensor. EF estimates were produced for the years 2010 to 2012 and locally validated with in-situ eddy-covariance data. The utility and robustness of the proposed indicator for operational monitoring of crop conditions at regional scale was assessed by comparing seasonal synthesis of EF and NDVI with official yield data.

Study Area Description
The study was focused on agricultural areas of Lombardy, in Northern Italy. The AOI is a part of Po valley and spans about 12000 km 2 . Topography of the study area is mostly flat, with altitude slowly degrading from about 150 m in the North West to about 20 m in the Eastern area. It is intensively cultivated with summer crops (i.e., corn (Zea mays L.), forage and rice (Oryza sativa L.)) and winter wheat (Triticum L.), which cover about 68% and 8% of the agricultural land, respectively. It is divided in 45 agronomical districts, showing marked differences in the main cultivated crops ( Figure 1).   Average air temperature is 13-14 • C, with summer (June-August) and winter (December-February) average values of 23-24 • C and 2-3 • C, respectively [51]. Precipitations are concentrated in autumn and spring, with a mean annual total of 800 mm. In this region, which is one of the most urbanized and industrialized of Italy, water management is critical during the dry summer months when multiple and conflicting usages (i.e., civil, industrial, agriculture and hydroelectric) can reduce water availability for irrigation. Official statistics of annual crop yield over the agricultural districts are publicly made available by the Lombardy Environmental Protection Agency (ARPA).
As the area shows flooded, irrigated and not-irrigated crops, is considered suitable for regional monitoring of crop/surface moisture conditions with the triangle approach, which requires a large number of pixels to cover a wide range of soil moisture and vegetation conditions [11,23]. In the last decade, the area experienced water scarcity in summer months when water availability could not satisfy the multiple water requirements [52] due to a lack of precipitations in key periods of crops growing cycle and/or reduced snowfall in winter. Irregular precipitations might also be accompanied by excessively high temperatures, such as in the year 2003, when corn production in Italy was 36% lower than the 1998-2002 average [53]. Besides climatic conditions, non-optimal timing in water supply to farmers by public authorities can worsen crop conditions, as in 2012 when corn production decreased by about 15-20% due to water shortage [54]. The reduction in corn production was worsened by the outbreak of aflatoxin contamination in corn stocks [55] induced by the occurrence of water stress conditions during critical phenological stages.
To face these problems, regional authorities involved in agro-monitoring over this area require therefore detailed information about the temporal and spatial variations of crop conditions. Although periodical agro-meteorological bulletins which analyze the impact of water stress on crop production are already issued for the study area by international agencies (e.g., Joint Research Centre-Monitoring Agricultural Resources Unit-JRC-MARS [56]; Global Agricultural Monitoring Initiative-EOGLAM [57]), they are not sufficiently detailed in terms of spatial resolution. More specific regional monitoring programs providing dedicated analyses for the agricultural areas of Lombardy are therefore needed.

Satellite Data and Preprocessing
MODIS 250 m resolution 16-day Vegetation Indices (VI) composites (Products MOD13Q1/ MYD13Q1) and 1 km Land Surface Temperature and Emissivity (LST/E) daily (Product MYD11A1) data and the corresponding Quality Assurance (QA) information [58] over tile h18v04 and the period 2010-2012 have been downloaded and pre-processed using the MODIStsp "R" package [59].
Concerning vegetation indexes data, the use of composite data allowed to exploit VI already derived from the best possible observations in the composite period, thus reducing oscillations related to varying observation geometries and atmospheric conditions. Moreover the synergic use of Terra and Aqua products VIs, with a theoretical 8 days shift, allowed to analyze a denser time series.
The raw M*D13Q1 time series were smoothed in order to remove residual noise in the NDVI signal, using a two-step approach based on (i) outlier removal and (ii) weighted Savitzky-Golay filtering with weights assigned on the basis of QA information, considering the real date of observation derived from the Day of Acquisition MODIS layer. In particular, weights are computed based on pixel reliability and VI usefulness indicators derived from the MODIS QA layers, complemented with specific thresholds blue reflectance [60,61]. The algorithm output consists in a smoothed 250 m NDVI time series with an 8-day regular time step.
This product was finally resampled to 1 km spatial resolution through average aggregation, and daily time series were generated through linear interpolation between values of the nearest dates. This allows to reduce potential biases during the periods of rapid crop growth.
The Aqua MYD11A1 LST/E product provided 1 km daily surface temperature (T s ) accompanied by pixel-level QA flags. QA parameters were summarized in a three-level categorical Quality Flag (QF): QF1 (best quality) if all QA parameters reported the highest quality; QF0.5 (intermediate quality) if at least one of the quality layers reported a sub optimal estimation; and QF0 (bad quality) for not computed T s . QF layer was first exploited to investigate the influence of QA level on the accuracy of MODIS T s and afterwards to discard low quality pixels [62]. No smoothing process was applied to T s time series since small daily fluctuations of surface temperature may be related to actual changes rather than noise and a noise-reduction technique might cause loss of valid data. Moreover, temporal composites of T s cannot be used because the triangle method relies on the assumption of uniform atmospheric conditions and incoming energy over the study area, which cannot be assured over longer time periods [23].

Air Temperature Data from Meteorological Stations
Seventy-five automatic meteorological stations belonging to the Lombardy region network (Rete Regionale di Rilevamento Meteorologico) and managed by the Lombardy Environmental Protection Agency (ARPA Lombardia) are located within the study area (purple dots in Figure 1). Measurements acquired by the meteorological network are publicly available through the ARPA web site (http: //goo.gl/zqkmf0-last access May 2017). In particular, air temperature (T a ) is available with half-hour frequency. Daily T a measurements simultaneous with Aqua satellite overpass (1.30 pm) for the 2010-2012 period were downloaded and spatially interpolated the study area by applying an Inverse Distance Weighted (IDW) approach to derive daily 1 km maps.

Eddy Covariance Measurements
Two Eddy Covariance (EC) towers are located within the study area, respectively in Landriano (45.32 N; 9.27 E, 89 m a.s.l.) and Livraga (45.19 N, 9.57 E, 61 m a.s.l.) ( Figure 1-green diamonds). The towers are located within corn fields of about 10 and 12 ha [63,64] (see Appendix A). They are equipped with a three-dimensional sonic anemometer to measure sonic temperature and the three components of wind speed, an open-path gas analyzer to measure water vapor and carbon dioxide concentrations, a net radiometer and a thermo-hygrometer. Soil moisture, soil temperature and soil heat flux observations are also measured at different depths (10-30-50 cm, 4-8 cm and 6 cm respectively) every 30 min.
Data collected for the years 2010-2011 at Landriano and 2010-2012 at Livraga were used as ground truth information. In particular, air and surface temperature measurements at the EC towers were used to assess the accuracy of the T a and T s datasets used as inputs to compute EF, while EF estimates simultaneous to the Aqua satellite overpass (1.30 pm) were compared to satellite-based EF estimates.
To compute EF at the towers location, raw EC data were corrected to account for instrumental (e.g., axis rotation for tilt corrections, spike removal, time lag compensation) and physical errors (spectral information losses, air density fluctuations and air humidity effects on sonic temperature) following procedures from the literature [65], as implemented in the PEC (Polimi Eddy Covariance) software [64]. After these corrections procedures, the energy budget has been analyzed, showing a closure ranging between 0.92 and 0.97. According to literature, a lack of energy balance closure is common in EC measurements [66,67], and can be considered acceptable up to values of 0.6 [68]. Following the procedure developed by Twine et al. [69], latent and sensible heat fluxes were analyzed respecting the Bowen-ratio method to reach the energy balance closure. Net radiation (R n ), latent (λE), sensible (H) and soil (G) heat fluxes were then derived from corrected EC measurements and exploited to compute EF with a semi-hourly time step using the formulation of the energy balance [70]:

Ancillary Datasets
The Carta Uso Agricolo Annuale (CUUA) is a 20 m resolution land-use map produced each year by the Regional Agency for Agriculture and Forest Services of Lombardy (ERSAF) (http: //goo.gl/GJPU3U, last access May 2017) by combining farmers annual declarations with the official regional land use map (Destinazione d'Uso dei Suoli Agricoli e Forestali-DUSAF). The map includes 21 agricultural land-cover/use classes of which the most common over the study area are shown in Figure 1. The CUUA map was used to create a 1 km mask of arable land, by retaining only pixels including at least 60% of crop land and an average elevation below 100 m a.s.l. (corresponding to maximum elevation of irrigation systems in Lombardy plains). The arable land mask was used to discard from the analysis non-agriculture (e.g., urban, forest) and mixed land-use pixels (corresponding to about 37% of the area), which inclusion would have affected EF retrieval performance [23].

Computation of EF from MODIS Data
Methods for computing EF from satellite data generally rely on the evidence that the combination of T s and NDVI is diagnostic of surface moisture condition [37]. In particular, in this study we used the approach proposed by Jiang and Islam [16], which is based on the relationship between NDVI and the air-surface temperature difference (∆T = T s − T a ). The ∆T vs. NDVI bidimensional space exhibits for a given day a triangular shape, bounded by the wet (maximum evapotranspiration) and dry (minimum evapotranspiration) edges [71] (Figure 2). The wet edge position corresponds to the x-axis (∆T = 0) as established by the formulation of the latent heat flux [16]. The dry edge position is instead estimated as the regression line of maximum observed ∆T at discrete intervals of NDVI (green dots in Figure 2).

Ancillary Datasets
The Carta Uso Agricolo Annuale (CUUA) is a 20 m resolution land-use map produced each year by the Regional Agency for Agriculture and Forest Services of Lombardy (ERSAF) (http://goo.gl/GJPU3U, last access May 2017) by combining farmers annual declarations with the official regional land use map (Destinazione d'Uso dei Suoli Agricoli e Forestali-DUSAF). The map includes 21 agricultural land-cover/use classes of which the most common over the study area are shown in Figure 1. The CUUA map was used to create a 1 km mask of arable land, by retaining only pixels including at least 60% of crop land and an average elevation below 100 m a.s.l. (corresponding to maximum elevation of irrigation systems in Lombardy plains). The arable land mask was used to discard from the analysis non-agriculture (e.g., urban, forest) and mixed land-use pixels (corresponding to about 37% of the area), which inclusion would have affected EF retrieval performance [23].

Computation of EF from MODIS Data
Methods for computing EF from satellite data generally rely on the evidence that the combination of Ts and NDVI is diagnostic of surface moisture condition [37]. In particular, in this study we used the approach proposed by Jiang and Islam [16], which is based on the relationship between NDVI and the air-surface temperature difference (ΔT = Ts − Ta). The ΔT vs. NDVI bidimensional space exhibits for a given day a triangular shape, bounded by the wet (maximum evapotranspiration) and dry (minimum evapotranspiration) edges [71] (Figure 2). The wet edge position corresponds to the x-axis (ΔT = 0) as established by the formulation of the latent heat flux [16]. The dry edge position is instead estimated as the regression line of maximum observed ΔT at discrete intervals of NDVI (green dots in Figure 2). Given the position of a generic point i in the scatterplot, EF is estimated on the basis of its distance from the dry and wet edges [14]: Given the position of a generic point i in the scatterplot, EF is estimated on the basis of its distance from the dry and wet edges [14]: where ∆T i is the difference between surface (T s ) and air (T a ) temperature for pixel i; ∆T H (NDV I i ) is the ∆T of the dry edge in correspondence of the NDVI value of pixel i. Therefore, the closer the pixel is to the dry edge, the lower is EF (down to zero), meaning that it is characterized by low moisture conditions. The described method was used to create daily 1 km EF maps (EF d ) over the study area for 2010-2012. Daily NDVI values, derived as described in (Section 2.2), and daily ∆T maps derived as the difference between MODIS T s maps (Section 2.2), and T a maps derived from weather stations (Section 2.3.1) were used to build the ∆T. vs. NDVI daily scatterplots, and successively derive daily EF values for each (cloud free) pixel. Computation was not conducted on dates with more than 40% of cloud cover since they wouldn't include the full range of variability of thermal and crop cover conditions that is one of the essential element to reliably compute the dry and wet edges and therefore to apply the triangle method. To avoid problems due to residual cloud contamination, pixels below the lower 15th percentile of the NDVI distribution were not considered in the dry edge computation. The number and positions of NDVI intervals used to compute the regression line were determined with a logarithmic function, as proposed by Sturges [72].
Equation (2) provides instantaneous EF estimation at the time of satellite overpass. Although the diurnal behavior of EF exhibits a typical light concave-up shape with a minimum around noon [73], its variations are much lower compared to latent and sensible heat fluxes under clear sky conditions [17,73,74]. EF is therefore frequently assumed to be reasonably constant during daytime, and satellite-based EF estimates are often considered representative of overall daily conditions [23]. Peng et al. [17] added that the assumption of daytime constant EF is more robust if the estimation is performed around midday (as the case of MODIS Aqua products) and in clear sky conditions only (as is the case after the removal of contaminated data thanks to QF). Therefore, given that the aim is to estimate a weekly indicator of surface moisture conditions rather than to quantitatively estimates daily water needs, satellite based instantaneous EF maps were considered a suitable indicator of surface moisture conditions over large areas.
EF d maps were averaged to derive a weekly indicator of crop water stress (EF w ). The use of weekly averages is more suitable for regional crop monitoring since it reduces spatiotemporal gaps due to cloud cover and it decreases the influence of day to day variability due to external factors, e.g., wind speed and other atmospheric interferences [75]. A weekly indicator is also better suited for providing information for agro-meteorological bulletins, which are usually issued with a monthly or biweekly frequency during the crop season.

Comparison with In Situ Data
Comparison with in-situ data collected at the EC towers was conducted to verify the accuracy of input temperature maps (T a and T s ), and EF w estimates over the 1 km pixels corresponding to the two EC towers. In the case of EF, the weekly average was computed starting from the daily half-hourly EC estimations obtained at the hour of satellite overpass (13:30).
Averaging over the week allows the reduction of the impact of missing measurements due to instrumental and physical errors in data acquisition. Linear regression parameters (slope, intercept and coefficient of determination r 2 ) and indices of agreement (Mean Error (ME), Root Mean Squared Error (RMSE) and Relative Root Mean Squared Error (rRMSE)) were considered as performance indicators.
The comparison was limited to the growing period of the more common summer crops, i.e., between sowing (mid-April) and senescence (mid-August). Outside this period the estimation of EF is not reliable because the full range of thermal and vegetation cover conditions needed to calculate the dry edge are not fulfilled [6,23,32,76]. Moreover, frequent cloud cover in winter and autumn may prevent the use of satellite data for operational monitoring. Moreover, in the study area, information on crop conditions are necessary for water management and irrigation planning only in the hot and dry summer months, when water availability is a key issue. On the contrary, the autumn/winter season is not crucial for operational crop monitoring.

Comparison to Crop Yield Statistics
The usefulness of EF w as a spatially distributed indicator of surface moisture conditions was assessed over five agronomical districts of the study area. The selected districts are representative of the major crop types: corn (C1: Codogno plain, C2: Central Bresciana plain, and C3: Cremona plain), rice (R1: Western Lomellina) and forage (F1: Western Oltrepo Mantovano). EF w profiles over these five districts were analyzed to investigate whether the indicator (i) depicts the expected trend of crop seasonal dynamics in relation to known crop characteristics and agro-practices, and (ii) provides evidence of inter-annual variability related to water stress.
Corn districts C2 and C3 witnessed a significant decrease in production in 2012, whereas corn district C1 did not suffer from water stress due to sufficient irrigation. Cumulated summer EF w (June to August) over these three districts was calculated to better investigate whether EF w estimates can provide insights on the occurrence of summer droughts and their impact on crop production. Cumulated NDVI over the same period was also computed since it is a proxy of vegetation biomass frequently used to identify crop production shortages [77]. One-way analysis of variance (ANOVA) followed by the post-hoc Tukey HSD (Honestly Significant Difference) [78] multiple range test were carried out over the three districts in order to evaluate the diagnostic capability of both indicators in detecting differences between the years. Inter-annual variability of the indicators was finally compared with official statistics of corn yield at district level. dry summer months, when water availability is a key issue. On the contrary, the autumn/winter season is not crucial for operational crop monitoring.

Comparison to Crop Yield Statistics
The usefulness of EFw as a spatially distributed indicator of surface moisture conditions was assessed over five agronomical districts of the study area. The selected districts are representative of the major crop types: corn (C1: Codogno plain, C2: Central Bresciana plain, and C3: Cremona plain), rice (R1: Western Lomellina) and forage (F1: Western Oltrepo Mantovano). EFw profiles over these five districts were analyzed to investigate whether the indicator (i) depicts the expected trend of crop seasonal dynamics in relation to known crop characteristics and agro-practices, and (ii) provides evidence of inter-annual variability related to water stress.
Corn districts C2 and C3 witnessed a significant decrease in production in 2012, whereas corn district C1 did not suffer from water stress due to sufficient irrigation. Cumulated summer EFw (June to August) over these three districts was calculated to better investigate whether EFw estimates can provide insights on the occurrence of summer droughts and their impact on crop production. Cumulated NDVI over the same period was also computed since it is a proxy of vegetation biomass frequently used to identify crop production shortages [77]. One-way analysis of variance (ANOVA) followed by the post-hoc Tukey HSD (Honestly Significant Difference) [78] multiple range test were carried out over the three districts in order to evaluate the diagnostic capability of both indicators in detecting differences between the years. Inter-annual variability of the indicators was finally compared with official statistics of corn yield at district level. The reliability of the exploited NDVI, showed by the smoothed lines in Figure 3, is significantly increased by applying the multi-temporal smoothing process described in Section 2.2; yet this kind of filter cannot be used to Ta and Ts time series since daily fluctuations may be related to changes in The reliability of the exploited NDVI, showed by the smoothed lines in Figure 3, is significantly increased by applying the multi-temporal smoothing process described in Section 2.2; yet this kind of filter cannot be used to T a and T s time series since daily fluctuations may be related to changes in surface/air temperature and a noise-reduction technique might cause a loss of valid information. Scatter plots between in situ and estimated air and surface daily temperature are shown in Figure 4 and regression parameters are summarized in Table 1. Globally, estimated T a shows a very good accordance with data measured at the EC towers (r 2 = 0.95; ME = 1.37) with the exception of the Landriano site in 2010 (red dots in Figure 4a), which shows a noticeable overestimation (ME = 2.94, Table A1 in Appendix B). This good agreement is favored by the spatial homogeneity of air temperature over flat areas such as the Lombardy plain [80]. Moreover, both towers are surrounded by a very dense meteorological network with at least five stations in a range of 15 km around each tower (purple dots in Figure 1) thus making T a interpolation robust. surface/air temperature and a noise-reduction technique might cause a loss of valid information. Scatter plots between in situ and estimated air and surface daily temperature are shown in Figure 4 and regression parameters are summarized in Table 1. Globally, estimated Ta shows a very good accordance with data measured at the EC towers (r 2 = 0.95; ME = 1.37) with the exception of the Landriano site in 2010 (red dots in Figure 4a), which shows a noticeable overestimation (ME = 2.94, Table A1 in Appendix B). This good agreement is favored by the spatial homogeneity of air temperature over flat areas such as the Lombardy plain [80]. Moreover, both towers are surrounded by a very dense meteorological network with at least five stations in a range of 15 km around each tower (purple dots in Figure 1) thus making Ta interpolation robust.   The comparison of satellite Ts vs. EC measurements showed less satisfactory results, highlighting a consistent underestimation. The analysis of the impact of the quality of satellite observations, as reported in the QA layer of the satellite datasets (Section 2.2), on the accuracy of estimated surface temperatures highlighted that data flagged as good quality (QF1) provide a significantly better estimation (ME = −0.48 °C, r 2 = 0.74) compared to QF0.5 data (ME = −2.38 °C, r 2 = 0.74). This testifies that QF0.5 data, although being considered as intermediate quality, are still strongly affected by residual cloud cover and/or atmospheric interference and led us to decide to retain only QF1 data for EFd computation (See Table A1 in Appendix B for a more detailed analysis on Ts data quality through sites and years).

Comparison of Air and Surface Temperature with In Situ Measurements
This decision caused a reduction of the cardinality of available Ts estimates: for example, on the two pixels corresponding to the towers, out of the 1825 available dates only 36% was assigned to the higher quality level (QF1). This percentage raises to about 53% on summer months, suggesting that a much higher reliability can be achieved in EFd estimation in that period (See Figure A3 and A4 in Appendix B for a more detailed analysis on temporal variability of Ts data quality).  The comparison of satellite T s vs. EC measurements showed less satisfactory results, highlighting a consistent underestimation. The analysis of the impact of the quality of satellite observations, as reported in the QA layer of the satellite datasets (Section 2.2), on the accuracy of estimated surface temperatures highlighted that data flagged as good quality (QF1) provide a significantly better estimation (ME = −0.48 • C, r 2 = 0.74) compared to QF0.5 data (ME = −2.38 • C, r 2 = 0.74). This testifies that QF0.5 data, although being considered as intermediate quality, are still strongly affected by residual cloud cover and/or atmospheric interference and led us to decide to retain only QF1 data for EF d computation (See Table A1 in Appendix B for a more detailed analysis on T s data quality through sites and years).
This decision caused a reduction of the cardinality of available T s estimates: for example, on the two pixels corresponding to the towers, out of the 1825 available dates only 36% was assigned to the higher quality level (QF1). This percentage raises to about 53% on summer months, suggesting that a much higher reliability can be achieved in EF d estimation in that period (See Figures A3 and A4 in Appendix B for a more detailed analysis on temporal variability of T s data quality).

Comparison of EF w with In Situ Data
Scatterplots comparing EF w estimates derived from satellite and from EC data are shown in Figure 5, grouped by site and year; agreement metrics are summarized in Table 2.

Comparison of EFw with In Situ Data
Scatterplots comparing EFw estimates derived from satellite and from EC data are shown in Figure 5, grouped by site and year; agreement metrics are summarized in Table 2.  Satellite EFw values are generally lower than EC derived EFw, with an overall ME of −0.15, best (worst) performance are observed over Landriano 2011 (Livraga 2010). These results are consistent with Lu et al. [81], but opposite in sign to those obtained by de Tomás et al. and Peng and Loew [23,82]. This uncertainty in EF estimation may be related to the difference between the 1 km pixel size of MODIS maps used as inputs and the footprint of the EC towers (120 m 2 in Landriano and 140 m 2 in Livraga) [64]. The area surrounding the towers is fragmented and heterogeneous in terms of crop type: on average, corn covers 36% and 21% of the 1 km MODIS pixels at Landriano and Livraga, respectively. In these conditions, the 1-km Ts and NDVI values derived from MODIS are influenced also by the characteristics of the areas surrounding the EC towers, thus making the comparison with field-estimated EF difficult. This is further exacerbated by the fact that the ratio between the areas cultivated with winter and summer crops is variable from year to year ( Figure A1 in Appendix A). For example, in the case of Livraga 2012, 250 m NDVI time series highlight the effect of neighboring winter wheat cultivations on the satellite signal (Figure 3-bottom panel). This is further evident on the 1 km NDVI time-series, which in 2012 are found to be quite different from what would be expected for a summer-crop area ( Figure A2 in Appendix A).  Satellite EF w values are generally lower than EC derived EF w , with an overall ME of −0.15, best (worst) performance are observed over Landriano 2011 (Livraga 2010). These results are consistent with Lu et al. [81], but opposite in sign to those obtained by de Tomás et al. and Peng and Loew [23,82]. This uncertainty in EF estimation may be related to the difference between the 1 km pixel size of MODIS maps used as inputs and the footprint of the EC towers (120 m 2 in Landriano and 140 m 2 in Livraga) [64]. The area surrounding the towers is fragmented and heterogeneous in terms of crop type: on average, corn covers 36% and 21% of the 1 km MODIS pixels at Landriano and Livraga, respectively. In these conditions, the 1-km T s and NDVI values derived from MODIS are influenced also by the characteristics of the areas surrounding the EC towers, thus making the comparison with field-estimated EF difficult. This is further exacerbated by the fact that the ratio between the areas cultivated with winter and summer crops is variable from year to year ( Figure A1 in Appendix A). For example, in the case of Livraga 2012, 250 m NDVI time series highlight the effect of neighboring winter wheat cultivations on the satellite signal ( Figure 3-bottom panel). This is further evident on the 1 km NDVI time-series, which in 2012 are found to be quite different from what would be expected for a summer-crop area ( Figure A2 in Appendix A).    Figure 6 shows examples of daily T a , T s , ∆T, NDVI and EF d maps for 10 representative clear-sky dates in 2010. No maps are showed for December and November because no cloud-free data were available. T a and T s maps highlight the climatic variability of the study area, with cold winter months, mild spring and hot summer; ∆T maps show the highest values between April and May, in relation to the presence of vast bare soil areas associated with increasing incoming radiation and temperature. NDVI maps depict the growing season of the major summer crops, with the peak of vegetation greenness between June and August (NDVI > 0.5). They also highlight areas where forage and winter wheat are cultivated (October to February), mainly in the central/eastern regions of the study area. The fifth column of Figure 6 shows the ∆T vs. NDVI scatterplots derived from the input maps and used for EF d computation, with the dry-edge shown as a green line. It can be observed that the dry-edge slope and position show a strong seasonal variability due to the succession of warm months with pixels characterized by ∆T above 10 • C and cold months when the difference between soil and air temperature is usually low (see Appendix C for a more detailed analyses of dry-edge seasonal variability).

Spatially Distributed Estimation of EF d
EF d maps are shown in the last column, with lower EF d areas (drier crop conditions) in orange/red colors, and high EF d area (well-watered areas) shown in blue colors. Higher EF d values are obtained in summer months, when evapotranspiration is maximum and the well-watered full vegetated pixels are close to the wet edge in the ∆T vs. NDVI scatter plots. Districts where rice is the main crop (see CUAA map in Figure 1) are constantly characterized by high EF d values in summer, because paddies are flooded and never face water stress. Similar summer patterns are obtained in most of the corn districts, because corn is usually well irrigated by the widespread irrigation system of the region.
White areas inside the study area in Figure 6 are due to masked pixels (arable land mask or cloudy conditions), which is more frequent in January to April. Low quality pixels can also be present, to a less extent, during the summer season (Appendix B) thus reducing the available number of reliable EF d estimations. The use of an aggregated weekly moisture indicator (EF w ), reduces the likelihood of missing information and is therefore more suited for analyzing the overall seasonal trend of crops conditions. This also allows the comparison between different years to highlight inter-annual variations due to unfavorable years characterized by water scarcity or biotic/abiotic stress. Figure 7 shows EF w and average weekly NDVI over five agronomical districts in the study area (bold red polygons in Figure 1) from 2010 to 2012.

Spatiotemporal Patterns of EF Weekly Indicator
NDVI and EF profiles are smoothed with a Generalized Additive Model (GAM) in order to better highlight seasonal trends. Both NDVI and EF w show clear seasonality in the districts dominated by rice (R1) and corn (C1, C2 and C3), with highest values between July and August. Only district F1 (Western Oltrepo Mantovano) shows no seasonal behavior due to the dominance of meadows and Alfalfa (Medicago sativa L.) fields (green areas in the CUUA map). These perennial cultivations have a long growing cycle and they are harvested for fodder more than once in a year [83]. In this area, the variability in forage management, the lack of a widespread irrigation system and the presence of winter wheat (red spots in the CUUA map) are the reasons for lower summer NDVI and EF compared to the other districts.
The rice district (R1) shows the highest summer values of NDVI and EF w , and the lowest inter-annual variability; indeed, rice does not face water stress because it is traditionally grown in flooded conditions [84]. It also shows lower values at the end of the season compared to corn districts, where after corn harvesting stubbles are left over the fields, and plants tend to re-germinate thus increasing NDVI and EF. Moreover, rice is typically cultivated as single crop [85] thus leading to lower values of NDVI and EF at the beginning of the season. EF w time series over the corn districts (C1, C2 and C3) show high values during the summer season since corn is usually well irrigated and hardly faces water stress. EF w is however lower and more variable from year to year than in the rice district. In particular, for C2 (Central Bresciana plain) and C3 (Cremona plain) districts EF w values after July are lower in 2012 compared to other years ( Figure 7). These two districts were characterized in 2012 by low and erratic rainfall that affected crop production. Similar results were achieved by Abbas et al. [86] in China by analyzing a multi-temporal spatially averaged drought indicator obtained from MODIS T s and NDVI.
White areas inside the study area in Figure 6 are due to masked pixels (arable land mask or cloudy conditions), which is more frequent in January to April. Low quality pixels can also be present, to a less extent, during the summer season (Appendix B) thus reducing the available number of reliable EFd estimations. The use of an aggregated weekly moisture indicator (EFw), reduces the likelihood of missing information and is therefore more suited for analyzing the overall seasonal trend of crops conditions. This also allows the comparison between different years to highlight interannual variations due to unfavorable years characterized by water scarcity or biotic/abiotic stress. Figure 7 shows EFw and average weekly NDVI over five agronomical districts in the study area (bold red polygons in Figure 1) from 2010 to 2012.   Cumulated EF shows areas with a major decrease from 2011 and 2010 to 2012, up to −30%. While the corn yield difference maps confirm this difference, no clear difference is appreciable in the NDVI ∆∑ Summer maps (Figure 8). This result suggests that summer cumulated NDVI is less sensitive to water stress periods compared to EF; indeed, only prolonged water stress conditions significantly affect crop biomass (i.e., NDVI). On the contrary, the EF indicator is more sensitive to transient unfavorable conditions since surface temperature promptly increases with respect to air temperature (due to rapid adjustments in stomata closure) if crop water content decreases due to decreased water availability [87].

Spatiotemporal Patterns of EF Weekly Indicator
Results of ANOVA and post-hoc Tukey test are showed in Table 3a. Average (µ) cumulated EF in 2012 over C2 and C3 corn districts is statistically lower compared to the other years (α = 0.001), suggesting that corn production has faced a remarkable yield loss due to local dry conditions and scarce water availability for irrigation. No decrease of cumulated EF is observed over the Codogno plain district (C1) in 2012. The post-hoc test does not highlight statistically significant differences of the cumulated NDVI in any of the districts. This suggests that EF might be a more sensitive indicator for identifying changes in crop production related to water shortages. Compared to 2001-2011, average 2012 corn yield shows a stable value over C1 and a significant change over C2 (−5 q/ha) and C3 (−21 q/ha) districts (Table 3b), in agreement with EF indicator.
Δ∑Summer maps (Figure 8). This result suggests that summer cumulated NDVI is less sensitive to water stress periods compared to EF; indeed, only prolonged water stress conditions significantly affect crop biomass (i.e., NDVI). On the contrary, the EF indicator is more sensitive to transient unfavorable conditions since surface temperature promptly increases with respect to air temperature (due to rapid adjustments in stomata closure) if crop water content decreases due to decreased water availability [87]. Results of ANOVA and post-hoc Tukey test are showed in Table 3a. Average (μ)cumulated EF in 2012 over C2 and C3 corn districts is statistically lower compared to the other years (α = 0.001),

Conclusions
In this work, the triangle method was applied to retrieve daily EF estimates and a weekly indicator of surface moisture status (EF w ) in Northern Italy. Input data include (i) daily T a maps derived from the interpolation of meteorological stations data, (ii) daily T s maps derived from MODIS Land-Surface Temperature product, and (iii) NDVI maps derived from temporal interpolation of 16-day MODIS Aqua and Terra vegetation indexes products. The main advantage of this method lies in its reliance on just few input spatial data, which makes it suitable for large areas and operative implementation.
Reliability of the estimations depends on AOI characteristics ( [23,43,44], see Section 1.1 for details) and accuracy of the input temperature data [21]. The analyzed area, thanks to its great variety of crops with different water managements allows the computation of reliable triangular shape scatterplots in summer, when drought events could occur.
Availability of T a data simultaneous to satellite overpass is fundamental to compute the wet edge. In our study, interpolated T a maps showed a more than satisfactory agreement with field measurements, thanks to the dense meteorological stations network available and by the flat topography of the area. These favorable conditions are not frequently met worldwide, but alternative sources for T a estimation could be exploited, such as using data provided by the European Centre for Medium-range Weather Forecasts (ECMWF) [2] or exploiting an empirical method based on T s , NDVI and albedo datasets as tested by Sun et al. [76].
The accuracy of MODIS-derived surface temperature (T s ) is greatest for high quality satellite observations (QF1) while intermediate-quality data provide underestimated T s . Discarding pixels with lower quality is therefore necessary to guarantee reliable EF estimates even if the removed data introduce spatial gaps. In many cases, data might not be enough to represent the variability of surface conditions required for estimating the dry edge's parameters, making EF estimation unreliable or impossible on some days. The use of a weekly indicator (EF w ) copes with these limitations, as averaging daily values reduces the impact of both noise in instantaneous estimates and gaps due to unfavorable atmospheric conditions. Moreover, the QA dataset showed that, in the dry spring/summer months, a sufficient number of good-quality observations was available for the computation of a reliable weekly-averaged EF. In particular, within the most crucial period for drought crises (from June to August), at least one clear sky satellite image every 3 days was on average available, with a peak of 90% of cloud free observations in July, therefore assuring operational conditions for crop monitoring.
Comparison of the EF w maps against in-situ data showed that estimates are generally in line with weekly-averaged EC measurements, with performance comparable with those of previously published works relying on coarse spatial resolution satellite data. Yet, the simplification of the triangle model, such as the assumption of homogeneous conditions of surface available energy at the time of satellite overpass, leads to a bias in EF estimates especially at the regional scale [21].
The agreement between MODIS-derived and field-measured EF is also dependent on the difference between MODIS spatial resolution and EC footprint: indeed, this difference could be enhanced in fragmented and heterogeneous agricultural landscapes, where several crop types combined with different agronomic practices are interspersed in a small territory. Moreover, the intra-pixel landscape heterogeneity can lead to EF uncertainties since the methodology involves an assumption of uniform aerodynamic properties and physiological behavior of the surface [23].
While more accurate estimates could be achieved by using higher spatial resolution thermal imagery (e.g., de Tomás et al. [23]), satellite sources for this kind of data are currently characterized by longer revisiting time (e.g., 16 days for Landsat 8), that is not sufficient for repeated monitoring along the crop season. Unfortunately, new satellite sensors such as the Sea and Land Surface Temperature Radiometer (SLSTR), onboard ESA Sentinel 3, has thermal infrared bands with the same spatial resolution as MODIS. This limitation however doesn't prevent the operational use of the EF w indicator, whose scope is to support management at regional/district scale with a spatially distributed proxy of surface moisture conditions, rather than to assess water needs at field level.
Indeed, at the agricultural district scale, EF w showed temporal variations consistent with the agro-practices for the major crops in the study area (rice, forages and corn), as well as being correlated to crop production variability induced by water scarcity (i.e., the reduction of corn yield in 2012). These analyses suggest that EF w can provide information on crop conditions complementary to NDVI, which is commonly used as a proxy of annual biomass production [88]. EF w could be a useful indicator for improving agro-meteorological bulletins providing information on crop status during the growing season at regional scale, allowing to early highlight unfavorable conditions for crop yield, and potentially to implement mitigation activities.
Future work will be dedicated to estimating and analyzing the proposed indicator over the entire MODIS archive (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) Table A1 shows indices of agreement between in situ and estimated air and surface daily temperature. Compared to Table 1, regression parameters are divided by site and year. Compared to Ts, Ta data show more homogeneous goodness of fit measures across sites and years. ME values show that Ts estimates labeled as good quality (QF1) are characterized by a lower underestimation error compared to low quality data (QF0.5) in all sites and years except Landriano 2010 (ME = 1.03 and 2.69 for QF0.5 and QF1 respectively). This station was the only case where eddy data were available in wintertime, when satellite estimation of Ts could be more problematic; yet, water stress monitoring is less meaningful in wintertime.   Table A1 shows indices of agreement between in situ and estimated air and surface daily temperature. Compared to Table 1, regression parameters are divided by site and year. Compared to Ts, Ta data show more homogeneous goodness of fit measures across sites and years. ME values show that Ts estimates labeled as good quality (QF1) are characterized by a lower underestimation error compared to low quality data (QF0.5) in all sites and years except Landriano 2010 (ME = 1.03 and 2.69 for QF0.5 and QF1 respectively). This station was the only case where eddy data were available in wintertime, when satellite estimation of Ts could be more problematic; yet, water stress monitoring is less meaningful in wintertime. Appendix B. Table A1 shows indices of agreement between in situ and estimated air and surface daily temperature. Compared to Table 1, regression parameters are divided by site and year. Compared to T s , T a data show more homogeneous goodness of fit measures across sites and years. ME values show that T s estimates labeled as good quality (QF1) are characterized by a lower underestimation error compared to low quality data (QF0.5) in all sites and years except Landriano 2010 (ME = 1.03 and 2.69 for QF0.5 and QF1 respectively). This station was the only case where eddy data were available in wintertime, when satellite estimation of T s could be more problematic; yet, water stress monitoring is less meaningful in wintertime. Table A1. Indices of agreement between estimated and measured T a and T s divided by sites and years; for surface temperature the accuracy metrics are provided for high (QF1) and low (QF0.5) quality levels.  Figure A3 shows the proportion of the daily quality flags falling in the three classes (see Section 2.2) for each month from 2010 to 2012 at the two EC towers. In winter months (November to February), more than 70% of daily T s data are missing (QF0) due to frequent cloud cover, but in the spring-summer period (March to September) at least 60% of the T s data per month are produced (QF0.5, QF1), with an even greater proportion in the driest months of the year, which are of key relevance for water stress monitoring (peak of 90% of cloud free dates in July). Cardinality of available T s data is consistent over sites and years.   Figure A3 shows the proportion of the daily quality flags falling in the three classes (see Section 2.2) for each month from 2010 to 2012 at the two EC towers. In winter months (November to February), more than 70% of daily Ts data are missing (QF0) due to frequent cloud cover, but in the spring-summer period (March to September) at least 60% of the Ts data per month are produced (QF0.5, QF1), with an even greater proportion in the driest months of the year, which are of key relevance for water stress monitoring (peak of 90% of cloud free dates in July). Cardinality of available Ts data is consistent over sites and years. Boxplot in Figure A4 shows the statistics of Ts satellite estimation per months for good (QF1) and low quality data (QF0.5). Either time series shows a coherent temporal behavior. However, the medians of the boxes show that QF0.5 data always have lower values because of the well-known atmospheric interference. The use of QF0.5 data could lead to a misleading computation of ∆T, as is the case. Boxplot in Figure A4 shows the statistics of T s satellite estimation per months for good (QF1) and low quality data (QF0.5). Either time series shows a coherent temporal behavior. However, the medians of the boxes show that QF0.5 data always have lower values because of the well-known atmospheric interference. The use of QF0.5 data could lead to a misleading computation of ∆T, as is the case.

Appendix C
Seasonal trends of the parameters of the regression line representing the dry edge are shown in Figure C1. Monthly statistics of the slope and intercept are plotted against time with a box-plot representation. The summer months, from June to August, when temperatures are higher and water shortage could became an issue, are highlighted by yellow boxes in the graphs and delimited by red dotted lines.
With the onset of the growing season in late spring and the beginning of the summer period, dry edge becomes gradually steeper reaching average minimum slope between July and August as also displayed by scatter plots in Figure 6. Dry edge lines with high intercept and low slopes indicates dates characterized by greater difference between air Ta and surface Ts temperature, a thermal condition typical of hot months when stressed pixels are present in the area together with wellwatered areas; in these months of the year a wider range of thermal conditions of crops is reached.
Winter and spring months show, on average, dry edges with smaller slopes and intercepts close to zero due to both lower absolute temperature values and smaller difference between surface and air temperature (ΔT). During this time of the year in the study area, atmospheric water demand, hence potential evapotranspiration, is very low, in fact winter crops are not irrigated. Clearly, months when crop biomass is low or vegetation absent, the full range of thermal conditions are not covered and the estimated dry edge could be far from the theoretical one [6,23,32,76].
(a) Figure A4. Values of estimated T s over the two eddy towers for months and QF.

Appendix C.
Seasonal trends of the parameters of the regression line representing the dry edge are shown in Figure A5. Monthly statistics of the slope and intercept are plotted against time with a box-plot representation. The summer months, from June to August, when temperatures are higher and water shortage could became an issue, are highlighted by yellow boxes in the graphs and delimited by red dotted lines.
With the onset of the growing season in late spring and the beginning of the summer period, dry edge becomes gradually steeper reaching average minimum slope between July and August as also displayed by scatter plots in Figure 6. Dry edge lines with high intercept and low slopes indicates dates characterized by greater difference between air T a and surface T s temperature, a thermal condition typical of hot months when stressed pixels are present in the area together with well-watered areas; in these months of the year a wider range of thermal conditions of crops is reached.
Winter and spring months show, on average, dry edges with smaller slopes and intercepts close to zero due to both lower absolute temperature values and smaller difference between surface and air temperature (∆T). During this time of the year in the study area, atmospheric water demand, hence potential evapotranspiration, is very low, in fact winter crops are not irrigated. Clearly, months when crop biomass is low or vegetation absent, the full range of thermal conditions are not covered and the estimated dry edge could be far from the theoretical one [6,23,32,76].
The year 2012 displays the summer months with the highest intercepts, in particular in July (median value for intercept equal to 21 • C) consistent with the particular dry weather conditions of this year. The seasonal behavior displayed in the boxplots fits with the typical average climate condition of the study area [2] and are consistent with results from previous work carried out in different ecoregions [20,[36][37][38]76,89,90].
Winter and spring months show, on average, dry edges with smaller slopes and intercepts close to zero due to both lower absolute temperature values and smaller difference between surface and air temperature (ΔT). During this time of the year in the study area, atmospheric water demand, hence potential evapotranspiration, is very low, in fact winter crops are not irrigated. Clearly, months when crop biomass is low or vegetation absent, the full range of thermal conditions are not covered and the estimated dry edge could be far from the theoretical one [6,23,32,76]. The year 2012 displays the summer months with the highest intercepts, in particular in July (median value for intercept equal to 21 °C) consistent with the particular dry weather conditions of this year. The seasonal behavior displayed in the boxplots fits with the typical average climate condition of the study area [2] and are consistent with results from previous work carried out in different ecoregions [20,[36][37][38]76,89,90].