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

Abstract

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 (EF_w) were then derived at regional scale for the years 2010, 2011 and 2012 as a proxy of overall surface moisture condition. EF_w 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 EF_w 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.

1. 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 distributed information at province/regional scale can hardly be extrapolated due to heterogeneity of land surface and complexity of the hydrological processes [3].

Since Remote Sensing (RS) data can provide a large variety of information on crop status and surface hydrological conditions, which are of key importance to highlight potential criticalities to support water planning and provide quantitative data for management [4,5], RS techniques are recognized as the only viable means to map ET at regional scale in a consistent and economically feasible way [6]. Instantaneous values of ET at satellite overpass can be used as diagnostics for surface status [7], or as controls for hydrological models through assimilation schemes [8]. Several approaches have been developed to estimate ET and/or indicators of water stress, spanning from simple empirical methods to more complex energy balance models, as fully described in literature reviews [9,10,11,12].

1.1. The Triangle Method: A Short Review on Past Applications and Recent Improvements

Among the available RS-based approaches for ET analysis, the triangle method, based on the work of Price in the early 1990s [13] and later elaborated by Jiang and Islam [14,15,16], has been widely exploited to estimate the Evaporative Fraction (EF), which is the ratio between the latent heat flux and the total available energy at canopy surface [5,17]. EF is both a key parameter to estimate ET as well as a direct indicator of surface moisture conditions [18] and water stress itself [19,20]. This approach, regarded as a simplification of more complex models such as the Surface Energy Balance Algorithm for Land (SEBAL) [21], is suitable for large area monitoring of surface moisture conditions [14,22]. Among its advantages are the simple parameterization/calibration and its reliance on operational satellite data [6,20,23], which allow to spatially explicitly estimate EF in near real-time over large areas.

In its original formulation, the triangle method builds on the triangular shape of the scatter plot of remotely sensed surface/canopy temperature (T_s) versus a Vegetation Index (VI) such as the Normalized Difference Vegetation Index (NDVI) [10,24]. This scatterplot is used to compute the so-called wet and dry edges: the wet edge corresponds to areas where EF is at maximum (i.e., maximum evapotranspiration) whereas the dry edge corresponds to areas where the EF is close to zero (i.e., limited water availability) [23]. Given two co-registered raster images containing temperature and VI value, EF can then be estimated for each pixel based on its relative position with respect to the two edges (see Section 3.1 for details).

Alternative formulations of the method were also proposed in later studies. In particular, using the difference between air surface and air temperature (ΔT = T_s − T_a) instead than simply T_s to build the aforementioned scatterplot was proposed by Moran et al. [25] and Jiang and Islam [16], and successively adopted by numerous studies (e.g., [26,27,28,29]). Additionally, other RS-derived indicators of vegetation characteristics (e.g., fractional cover, albedo, etc.) were proposed as an alternative to NDVI for computing the dry and wet edges [6,30].

The triangle method has been exploited starting from data acquired by several satellite platforms/sensors, among which NOAA-Advanced Very High Resolution Radiometer (AVHRR) [31], Meteosat Second Generation—Spinning Enhanced Visible and Infrared Imager (MSG SEVIRI) [32], NASA Moderate Resolution Imaging Spectroradiometer (MODIS) [2], Envisat-Advanced Along-Track Scanning Radiometer (AATSR) and Medium Resolution Imaging Spectrometer (MERIS) [23], Landsat Thematic Mapper/Enhanced Thematic Mapper (TM/ETM+) [23] as well as airborne sensors [28]. One of its main advantages is its suitability for hydrological studies where field data availability is limited or missing. Thanks to this, it was implemented in operational monitoring systems such as the EVapotranspiration Assessment from SPAce (EVASPA) tool [33], exploited to derive the Temperature Vegetation Dryness Index (TVDI) [25] used for vegetation water stress detection [34,35,36,37] and soil moisture estimation [38,39] in semi-arid, temperate and tropical areas [18] and used to improve the performances of hydrological models [40].

Nevertheless, it is important to remind that several issues should be addressed for a full exploitation of satellite based EF/ET estimations.

A first assumption of the method is that incoming energy, aerodynamic properties and atmospheric conditions have to be reasonably uniform over the area at the time of satellite overpass [2,23,32] and the AOI should be as flat as possible in order to identify a proper triangular shape in the pixel distribution [24]. Tave et al. [41] recently tried to solve this requirement by improving the triangle method to compute EF for different elevation zones rather than keeping the estimation of flat area only.

Second, a key issue for applying the triangle method is the proper computation of the wet and dry edges, which influences the accuracy of EF estimation [6]. The method can therefore be applied only under specific conditions, in relation to study area characteristics and EO data exploited, and a proper data handling is needed to retrieve reliable estimates. In particular: (i) a large range of combined soil moisture status and vegetation characteristics must be present in the study area to represent a wide range of evapotranspiration conditions [37,42], and (ii) this variability must be properly captured by the RS data used. These conditions influence the required extension of the area of interest (AOI) analyzed, which should be sufficiently large [43] to allow the construction of reliable scatterplots as a function of the spatial resolution of the EO data used. According to previous studies, the triangular shape depends more on the number of pixels, rather than on the spatial resolution [24,44]. Hence, for small AOIs (e.g., ~4000 km²) HR satellite data may allow obtaining more accurate EF estimations. On the other hand, while working with satellite data at coarse resolution the AOI should be bigger since a larger area is needed for an appropriate definition of the triangle shape [23].

If these conditions are not completely fulfilled, the shape of the temperature vs. VI scatterplot may be “flawed” (i.e., not triangular at all), leading the estimated dry and wet edges to be far from the theoretical ones, and therefore to large inaccuracies in EF estimation [32].

To deal with the aforementioned problem, many authors proposed modifications and advancements to the “standard” method allowing a better and more stable computation of the edges. Tang et al. [6] proposed a specific pre-processing technique for identifying outliers in the temperature vs. VI scatterplot space and removing them before the computation of the edges. De Tomas et al. [23] demonstrated how this approach can improve EF estimates when HR data (e.g., Landsat like) are used. Maltese et al. [28] compared scatterplots derived from different dates suggesting that, to cope with the uncertainness in the dry and wet edges identification, a multi-temporal analysis should be exploited to include outer extremes in soil water content. Minacapilli et al. [2] further developed this idea building the temperature vs. VI scatterplot in the temporal domains for each pixel of the AOI using all available dates; application of this approach resulted in an accurate estimation of ET at regional scale. Other recent contributions determined the theoretical edges on the base of the surface energy balance principle rather than by identifying the boundaries empirically in the data space; EF is then estimated for each pixel using the observed temperature vs. VI scatterplot [22,32]. Comparison with in situ data of soil moisture and EF revealed that this approach could increase the accuracy of the triangle method. However, its implementation requires a more complex parameterization and requires input field data not be always available over the AOI.

Another important issue discussed in literature is whether the output of the triangle method can be considered representative of the daily condition or the instantaneous momentum. Most approaches based on satellite data assumed EF to be relatively constant in daytime (“self-preservation of EF” [23]), so that satellite-based instantaneous EF estimates can be representative of daily conditions [30] (see Section 3.1 for details).

Nevertheless, for some scientific applications it is necessary to more rigorously upscale the instantaneous estimation to daily conditions taking into account the diurnal cycle using exogenous information such as meteo data [45]. To perform this daily upscaling, Trezza [46] proposed to multiply the instantaneous EF by the daily reference ET computed with FAO approach [47] while Ryu et al. [48] exploited the variations of the daily extraterrestrial solar radiation. Finally, Tang et al. [49] recently presented a promising approach for upscaling, based on half-hourly meteo data (i.e., air temperature, 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.

1.2. 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.

2. Materials

2.1. 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². 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.

2.2. 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].

2.3. Meteorological Data

2.3.1. 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.

2.3.2. 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]:

$$EF=\frac{\lambda E}{\lambda E+H}=\frac{\lambda E}{R_n-G}$$

(1)

2.4. 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].

3. Methods

3.1. 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).

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]:

$$E{F}_i=\left(\mathrm{\Delta }{T}_H\left(NDV{I}_i\right)-\mathrm{\Delta }{T}_i\right)/\left(\mathrm{\Delta }{T}_H\left(NDV{I}_i\right)\right)$$

(2)

where ΔT_i is the difference between surface (T_s) and air (T_a) temperature for pixel i; $\mathrm{\Delta }{T}_H\left(NDV{I}_i\right)$ 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 ($E{F}_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.

3.2. Accuracy and Usefulness Analysis

3.2.1. 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²) 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.

3.2.2. 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.

4. Results and Discussion

4.1. Comparison of Air and Surface Temperature with In Situ Measurements

Figure 3 shows a comparison between time-series of T_s and T_a estimated as described in Section 2.2 and Section 2.3.1 at the EC towers of Landriano and Livraga, and the corresponding measured values. Maximum T_s and T_a (about 30 and 35 °C, respectively) occur in the middle of July. This corresponds to the end of the vegetative period [79], as highlighted by the peak of the 250 m resolution NDVI time series (Figure 3—bottom panel).

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. Indices of agreement between estimated and measured T_a and T_s; for surface temperature the accuracy metrics are provided for high (QF1) and low (QF0.5) quality levels.

Variable	QF	N	Intercept	Slope	p value	r2	ME	RMSE	rRMSE
Ta [°C]	/	713	−1.36	1.00	<0.001	0.95	1.37	2.21	35.14
Ts [°C]	0.5	243	0.84	1.06	<0.001	0.74	−2.38	4.98	22.68
Ts [°C]	1	649	0.95	0.98	<0.001	0.74	−0.48	3.72	14.32

. Globally, estimated T_a shows a very good accordance with data measured at the EC towers (r² = 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. 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.

Data	QF	Site	year	N	Intercept	Slope	p Value	r2	ME	RMSE	rRMSE
Ta [°C]	/	Landriano	2010	233	−0.77	0.88	<0.001	0.96	2.94	3.61	17.33
			2011	109	−0.36	0.96	<0.001	0.99	1.35	1.42	31.92
		Livraga	2010	108	1.37	0.96	<0.001	0.98	−0.26	0.62	46.99
			2011	146	0.56	0.96	<0.001	0.99	0.53	0.71	47.64
			2012	117	0.36	0.95	<0.001	0.99	0.79	0.96	34.07
Ts [°C]	QF0.5	Landriano	2010	53	0.77	0.92	<0.001	0.85	1.03	3.50	15.68
			2011	60	7.88	0.88	<0.001	0.75	−4.73	5.32	23.60
		Livraga	2010	48	−1.63	1.14	<0.001	0.66	−2.54	4.36	18.75
			2011	46	8.47	0.86	<0.001	0.59	−4.59	5.57	24.53
			2012	36	5.95	0.96	<0.001	0.46	−4.80	7.62	43.46
Ts [°C]	QF1	Landriano	2010	112	0.42	0.88	<0.001	0.92	2.69	3.91	15.59
			2011	138	2.72	0.95	<0.001	0.67	−1.02	2.66	10.04
		Livraga	2010	114	3.37	0.92	<0.001	0.51	−0.87	3.56	14.13
			2011	147	7.99	0.84	<0.001	0.75	−3.05	3.77	14.15
			2012	138	6.51	0.78	<0.001	0.45	0.44	4.44	17.05

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.

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² = 0.74) compared to QF0.5 data (ME = −2.38 °C, r² = 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 Figure A3 and Figure A4 in Appendix B for a more detailed analysis on temporal variability of T_s data quality).

4.2. 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. Indices of agreement between weekly in situ and satellite EF_w. In bold the best values.

Station	Year	N	Intercept	Slope	r2	ME	RMSE
Landriano	2010	10	0.37	0.18	0.07	−0.18	0.22
Landriano	2011	14	0.22	0.56	0.41	−0.05	0.16
Livraga	2010	9	−0.32	0.79	0.52	−0.47	0.48
Livraga	2011	15	0.05	0.56	0.49	−0.21	0.28
Livraga	2012	10	0.23	0.44	0.18	−0.18	0.23
All	All	62	0.16	0.46	0.21	−0.2	0.28

.

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² in Landriano and 140 m² 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).

4.3. Spatially Distributed Estimation of EF_d

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).

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.

4.4. Spatiotemporal Patterns of EF Weekly Indicator

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.

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.

Figure 8 shows inter-annual difference, 2012 vs. 2010 and 2012 vs. 2011, of weekly NDVI and EF values cumulated from June to August (Δ∑_Summer) and maize yield for the corn agronomical districts. 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].

Results of ANOVA and post-hoc Tukey test are showed in

Table 3. (a) Results of ANOVA and post-hoc Tukey test (α = 0.001) for summer cumulated EF and NDVI among years (2010–2011) (n = cardinality; μ: district average value); (b) official annual corn yield statistics for the 2010–2012 period and last ten years average (2001–2011).

			Codogno Plain			Central Bresciana Plain			Cremona Plain
			C1			C2			C3
			2010	2011	2012	2010	2011	2012	2010	2011	2012
(a)	${{\displaystyle \sum }}_{Jun}^{Aug}E{F}_w$ [-]	Post-Hoc group	b	a	a	a	a	b	b	a	b
		μ	6.59	7.75	7.72	8.72	9	8.14	8.29	8.9	8.12
		n	229	232	228	283	359	355	231	255	248
	${{\displaystyle \sum }}_{Jun}^{Aug}NDV{I}_w$ [-]	Post-Hoc group	a	a	a	a	a	a	a	a	a
		μ	8.04	7.93	7.71	8.16	7.86	7.54	7.74	7.94	7.55
		n	231	240	228	131	263	230	235	264	248
(b)	Yield [q/ha]	Average 2001–2011		123			120			117
		Yearly	120	120	120	120	130	115	120	120	96

a. 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.

5. 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–2016) in order to assess its robustness for identification of other drought events occurred in the past decade, among which the ones of 2003 and 2015.