Modeling of Durum Wheat Yield Based on Sentinel-2 Imagery

: In this study, a modelling approach for the estimation/prediction of wheat yield based on Sentinel-2 data is presented. Model development was accomplished through a two-step process: ﬁrstly, the capacity of Sentinel-2 vegetation indices (VIs) to follow plant ecophysiological parameters was established through measurements in a pilot ﬁeld and secondly, the results of the ﬁrst step were extended/evaluated in 31 ﬁelds, during two growing periods, to increase the applicability range and robustness of the models. Modelling results were examined against yield data collected by a combine harvester equipped with a yield-monitoring system. Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) were examined as plant signals and combined with Normalized Difference Water Index (NDWI) and/or Normalized Multiband Drought Index (NMDI) during the growth period or before sowing, as water and soil signals, respectively. The best performing model involved the EVI integral for the 20 April–31 May period as a plant signal and NMDI on 29 April and before sowing as water and soil signals, respectively (R 2 = 0.629, RMSE = 538). However, model versions with a single date and maximum seasonal VIs values as a plant signal, performed almost equally well. Since the maximum seasonal VIs values occurred during the last ten days of April, these model versions are suitable for yield prediction.


Introduction
Durum wheat (Triticum turgidum L. var. durum) is a crop used for a variety of food products, mainly pasta. Almost 60% of the world durum wheat cultivated area, about 7.8 million hectares, is located at the Mediterranean basin [1]. Nevertheless, pasta consumption in the Mediterranean countries is higher than the local production, so food industries depend on imports from other durum-wheat-producing territories, mainly North America. Furthermore, the final production relies mainly on weather conditions, especially during the grain-filling stage [2]. The variability of the Mediterranean climate, exacerbated by the on-going climate change, causes great year-to-year fluctuations in durum wheat yields. This fact causes risks and uncertainty in the industry, grain marketing agencies, policymakers, and other involved entities, concerning the planning of their exports and imports. Early-season prediction of durum wheat yields is of vital importance for assisting the whole food production chain. Farmers as well, can adjust the farm inputs, such as fertilizers and irrigation, to meet the site-specific needs of the crop by implementing precision agriculture techniques, while the harvesting sector can plan its logistics by managing the harvester fleet and anticipating transport and storage requirements.
So far, yield prediction for winter wheat relies on either estimates and information gathered from experts or on outputs from a variety of crop simulation models, such as CERES [3], WOFOST [4], CROPSYST [5], or SAFY [6]. Even though models like the CERESwheat have been used quite successfully for more than 30 years [7], the main drawback is The criteria for selecting the fields were (i) the availability of yield map data and (ii) an adequate size, scheme, and orientation capable of providing a sufficient number of Sentinel-2, 10 × 10 m resolution pixels. All the fields were cultivated with durum wheat of different varieties (Iridae, Meridiano, Normano, Simeto, and Svevo) and belonged to numerous farmers who followed their own cultivation practices. That way, the model could be evaluated in different conditions and also, it could be ensured that it applied regardless of the cropping practices.    Raw yield data as provided by the harvesting machine (a) and after noise removal, rasterization, resampling, and collocation with Sentinel-2 pixels (b); Sentinel-2 NDVI on 29 April 2019 (c); soil ECa map spatial resolution after, rasterization, resampling, and collocation with Sentinel-2 pixels (d).
The 4.5 ha pilot field was cultivated with Svevo variety at a seed rate of 230 kg ha −1 . Six 30 × 30 m plots were accurately determined by GPS (Garmin, eTrex, Schaffhausen, Switzerland) to collocate with Sentinel-2 pixels (9 Senitnel-2 pixels per plot). All ecophysiological measurements were performed in these plots and compared with the corresponding Sentinel-2 data.

Field Measurements 2.2.1. Leaf Area Index
Leaf area index was measured with an AccuPAR LP-80 PAR/LAI Ceptometer (Decagon Devices, Inc., Pullman, WA, USA) following Norman and Jarvis model [38] of radiation transmission and scattering. Three LAI measurements were performed per plot, with 10 below and one above canopy PAR measurements per LAI measurement. PAR sampling points were randomly selected at 3 m intervals across three transects in each plot. For the below canopy measurements, the Ceptometer was placed at an angle of approximately 45 • in relation to cultivation lines. All measurements were performed during clear-sky days and around solar noon.

Water Potential
Water potential (Ψ) was measured using a Scholander-type pressure chamber (SKPM 1400, Skye Instruments Ltd., Liandrindrod, UK). Randomly selected whole plants across plots were wrapped in aluminum foil and sealed in plastic bags for 10 min and then cut and measured immediately with the pressure chamber.

Soil Electrical Conductivity
The soil apparent electrical conductivity (ECa) was measured using an EM38 instrument (Geonics Ltd., Mississauga, ON, Canada). The whole field was scanned at 22 March 2018 by walking the field in parallel lines, about 15 m apart, with the instrument handheld Agronomy 2021, 11, 1486 5 of 18 5 cm above the ground. The instrument records ECa at time intervals of 1s providing geotagged data used to create a point vector map. The vector map was initially interpolated into a raster with a resolution of 1x1m using the ordinary kriging process in the SAGA v. 7.9.0 free open-source GIS (SAGA User Group Association, 2020, Kent, UK). Accordingly, the raster map was resampled and collocated to the 10 × 10 m pixels of Sentinel-2. using the QGIS v. 3.16 open-source Geographic Information System.

Canopy Reflectance
Canopy reflectance was measured in situ with a portable JAZ spectroradiometer (OceanOptics, Dunedin, FL, USA), with spectral range from 350 to 1000 nm and spectral resolution approximately 0.3 nm. Measurements were taken with a 5 m length optical fiber in combination with a home-made supporting system at a height approximately 1 m from top of the canopy. The spectroradiometer was calibrated against a Spectralon ® white reflectance panel (Labsphere, Inc., North Sutton, NH, USA) before the measurements' onset and every 15 min, to track solar angle changes. All measurements were performed during clear-sky days and around solar noon time. At each plot 30 measurements were taken at 3-4 m intervals, by walking at three random lines forming a Z scheme.

Yield Measurement
At the end of the growing periods (June) the fields were harvested with a John Deere S660i combine harvester equipped with a yield mapping system providing, through the associated MyJD software, a yield map in a point vector format at a spatial resolution of 1.75 by approximately 2.5 m (depending on the travelling speed) (Figure 2a). The initial maps were further processed with the open-source Yield Editor software (v.2.0.7) [39] for the removal of outliers (due to start point, end point grain flow delays, swath width overlaps, rapid speed changes, etc.). Accordingly, the yield maps were interpolated to rasters by the IDW process of QGIS and finally, resampled at 10 × 10 m pixel size (Figure 2b) corresponding to the Sentinel-2 image pixels ( Figure 2c).

Satellite Data
A total of 61 cloud-free Sentinel-2 (A and B) images from October 2017 to June 2018 and October 2018 to June 2019 were downloaded from ESA's Copernicus Open Access Hub [40] (https://scihub.copernicus.eu/). The MultiSpectral Instruments (MSI), onboard the Sentinel-2 satellites, provide information at 13 spectral bands (443-2190 nm), at a variable spatial resolution from 10, 20, or 60 m pixel size and with a 5-day revisit time. In the present study, Level 2A (radiometrically and atmospherically corrected) bottom of atmosphere (BOA) reflectance products provided by ESA, were used. All images were resampled at 10 m pixel size using the SNAP-ESA Sentinels Application Platform v.6.0.4 [41] (http://step.esa.int) free open-source software. After data extraction for the pixels of the studied sites (221 pixels for the pilot field used for model development and 7436 pixels for the evaluation fields), time-series were constructed for the following vegetation indices: Normalized multiband drought index, NMDI = where R x , reflectance at wavelength x, with x denoting the center wavelength of the corresponding Sentinel-2 band.
Accordingly, the 221 + 7436 pixels × 4 VIs time series were linearly interpolated producing complete daily datasets for the corresponding growing periods (October to June). Linear interpolation and the absence of any smoothing process on the VI time series were chosen to keep processing as simple as possible, for the whole methodology to be easily applied in the modelling process. The interpolated daily datasets were subsequently used for the calculation of VI integrals for varying time intervals across the growing periods. Additionally, several soil indices were calculated for a single date before sowing (bare ground) and their formulations are shown in Table 1.

Statistics
The relationships between VIs and field-measured parameters were examined through correlation analyses, while the performance of estimation models was evaluated using linear and multilinear regression analyses, from which coefficient of determination (R 2 ) and root-mean-square error (RMSE) were derived (p-value < 0.05). All analyses were performed using JASP software v.0.14 [42].

Sentinel VIs and Field Measurements
Even though it is generally established that satellite VIs are good estimators of plant ecophysiological parameters, species-and site-specific relationships have to be well documented. To that purpose several Sentinel-2 VIs were examined against field-measured ecophysiological parameters in the pilot field.

Sensor Intercomparison
Canopy reflectance was measured with two different remote sensors-the MultiSpectral Instrument (MSI) onboard the Sentinel-2 satellite and the hyperspectral JAZ canopy spectroradiometer. Even though these two sensors have differences in spectral, spatial, and radiometric characteristics, their measurements may be used for the calculation of several widely used vegetation indices. In this study the JAZ spectroradiometer was used as a reference, since it has a super-fine spectral analysis (0.3 nm), it is calibrated before each measurement against a standard reflectance panel, and it is operated very close to the canopy avoiding any atmospheric disturbance effects. Accordingly, Sentinel-2 performance was evaluated on the basis of vegetation indices intercomparisons. For that purpose, the detailed spectral data of the JAZ spectroradiometer were averaged over the spectral range of the Sentinel-2 bands used in vegetation indices calculations. For each 30 × 30 m plot the JAZ measurements were averaged and compared to the average values of the corresponding nine Sentinel-2 pixels. Vegetation indices NDVI and EVI from Sentinel-2 correlated very well with JAZ indices, showing high correlation coefficients, as well as slopes and intercepts close to 1 and 0, respectively ( Figure 3). The large error bars for JAZ data may be considered rather reasonable due to its small measurement area (approximately 15 × 15 cm with the optical fiber at 1 m height above the canopy) resulting in detailed monitoring of the occurring variability, which is smoothed in the 10 × 10 m Sentinel-2 data.
Agronomy 2021, 11, x FOR PEER REVIEW 8 of 21 Figure 3. Relationships between JAZ and Sentinel-2 NDVI (a) and EVI (b). Data concern dates with both Sentinel-2 acquisitions and JAZ measurements. Each point corresponds to the mean from 10 JAZ measurements performed in one plot at one date and the mean from the corresponding 9 Sentinel-2 pixels ± standard deviation.

Leaf Area Index and Satellite VIs
Leaf area index was measured in the pilot field from February to May 2019. As shown in Figure 4, the seasonal course of the crop is well depicted in LAI, with low values during winter, gradually increasing until the April flowering period, followed by a rather steep decrease during the May seed-ripening period. This pattern was also followed by the satellite VIs ( Figure 5), resulting in very good correlations between field-measured LAI and VIs ( Figure 6).  Relationships between JAZ and Sentinel-2 NDVI (a) and EVI (b). Data concern dates with both Sentinel-2 acquisitions and JAZ measurements. Each point corresponds to the mean from 10 JAZ measurements performed in one plot at one date and the mean from the corresponding 9 Sentinel-2 pixels ± standard deviation.

Leaf Area Index and Satellite VIs
Leaf area index was measured in the pilot field from February to May 2019. As shown in Figure 4, the seasonal course of the crop is well depicted in LAI, with low values during winter, gradually increasing until the April flowering period, followed by a rather steep decrease during the May seed-ripening period. This pattern was also followed by the satellite VIs ( Figure 5), resulting in very good correlations between field-measured LAI and VIs ( Figure 6).
Agronomy 2021, 11, x FOR PEER REVIEW 2 of 21 drawback is that they require numerous data for weather, climate, soil, genotypes, and management practices that are often difficult to obtain. Moreover, these models are poor predictors when within-field spatial variability is considerable [8]. More recently, remotely sensed data have been used complementarily to improve the accuracy of crop models [9][10][11][12]. Remotely sensed data retrieved from satellite observations can provide quantitative and temporal information over large areas, through the use of vegetation indices (VIs), such as the normalized difference vegetation index (NDVI) [13], the enhanced vegetation index (EVI) [14], and the normalized difference red edge index (NDRE) [15]. Such VIs are closely related to vegetation biophysical parameters, such as the canopy leaf area index (LAI) or the fraction of absorbed photosynthetically active radiation (fAPAR) [16,17], which are commonly used in vegetation productivity modeling.
Even though spatial variability in crop productivity is highly associated with soil properties, water, and nutrient availability [2,18], remote-sensed data from the crop canopy have been proved a valuable tool for yield predictions, as they directly depict most of these effects on crop development [19]. Among the different approaches used for yield prediction, the empirical methods based on the straightforward application of a statistical regression between a VI and yield are the most common, because of their simplicity and limited data requirements [20][21][22][23][24][25][26]. The main drawback of this approach is that the relationships between VIs and yield are often limited to the regions for which they were calibrated and are not globally applicable [24,27]. Bhattacharya et al. [28] attributed the sitedependent nature of the spectral yield models on the saturation of some VIs (e.g., NDVI)  . Relationships between field-measured LAI and Sentinel-2 NDVI (a) and EVI (b). Data concern dates with both. acquisitions and field measurements. Each point corresponds to the mean from three LAI measurements performed in one plot at one date and the mean from the corresponding nine Sentinel-2 pixels ± standard deviation.   . Relationships between field-measured LAI and Sentinel-2 NDVI (a) and EVI (b). Data concern dates with both. acquisitions and field measurements. Each point corresponds to the mean from three LAI measurements performed in one plot at one date and the mean from the corresponding nine Sentinel-2 pixels ± standard deviation. Figure 6. Relationships between field-measured LAI and Sentinel-2 NDVI (a) and EVI (b). Data concern dates with both. acquisitions and field measurements. Each point corresponds to the mean from three LAI measurements performed in one plot at one date and the mean from the corresponding nine Sentinel-2 pixels ± standard deviation.

Plant Signal
Vegetation indices incorporating the red part of the spectrum (NDVI and EVI), where chlorophyll absorbs, were used as plant signal. Their performance was examined in the pilot field against final yield maps, resampled, and collocated according to Sentinel-2 spatial resolution (10 × 10 m), to enable a pixel-to-pixel comparison basis ( Figure 2). To that purpose, VI integrals of several time spans as well as single date VIs were compared to final yield data. Our initial hypothesis was that the whole growing period integral might show lower correlation with yield data compared to the integral for the period from flowering to harvesting, due to the fact that most of the assimilates are directed to seed formation after flowering [43]. Indeed, as shown in Table 2 and Figure 7, there is no correlation between the whole growing period integral (1 December-31 May) and yield (Figure 7a,b) for both NDVI and EVI. On the contrary, strong correlations appear when the integral from flowering to the end of the growing period (20 April-31 May, Figure 7c,d) or the best 20 days integral after the flowering period (6-25 May) are examined. Additionally, especially for NDVI, high correlations also appear when the single date with maximum seasonal value (max NDVI) is considered ( Table 2). In all high correlation cases NDVI seems to perform better than EVI. According to the above, only the high correlation cases from each VI are further considered in model configuration.  Water potential was measured in the pilot field only during April and May, s earlier dates plant shoots were too fragile to be inserted intact in the pressure bomb ratus. As shown in Figure 8, plants showed a good water status during April and ally lost water during the high temperature/low precipitation May. Noteworthy large error bars during the last May measurement, indicating a large in-field varia Figure 7. Relationships between final yield and VIs for the worst (a,b) and the best (c,d) integration periods. Data concern pilot field, corresponding to 221 pixels.

Water Signal
Water potential was measured in the pilot field only during April and May, since at earlier dates plant shoots were too fragile to be inserted intact in the pressure bomb apparatus. As shown in Figure 8, plants showed a good water status during April and gradually lost water during the high temperature/low precipitation May. Noteworthy are the large error bars during the last May measurement, indicating a large in-field variability. Subsequently, two satellite indices which incorporate water bands (NDWI and NMDI), were examined for their efficiency in depicting the water potential fluctuations. As shown in Figure 9, both indices showed good correlation with water potential. Even though NMDI is generally considered an advantageous water index as it incorporates two water bands, it showed weaker correlation with Ψ compared to NDWI. Subsequently, two satellite indices which incorporate water bands (NDWI and NMDI), were examined for their efficiency in depicting the water potential fluctuations. As shown in Figure 9, both indices showed good correlation with water potential. Even though NMDI is generally considered an advantageous water index as it incorporates two water bands, it showed weaker correlation with Ψ compared to NDWI. Subsequently, two satellite indices which incorporate water bands (NDWI and NMDI), were examined for their efficiency in depicting the water potential fluctuations. As shown in Figure 9, both indices showed good correlation with water potential. Even though NMDI is generally considered an advantageous water index as it incorporates two water bands, it showed weaker correlation with Ψ compared to NDWI. Figure 9. Relationships between water potential (Ψ) and Sentinel-2 water indices NDWI (a) and NMDI (b). Data concern dates with both Sentinel-2 acquisitions and field measurements. Each point corresponds to one Ψ measurement performed in one plot at one date and the mean from the corresponding 9 Sentinel-2 pixels ± standard deviation. Accordingly, in order to examine the importance of incorporation of a satellite water signal in the model, the correlation of the two water indices against final yield throughout the growing period were investigated. For both indices the most significant water signal (highest seasonal correlation) was found for 29 April (R = 0.641 and 0.550 for NDWI and NMDI, respectively), shortly after flowering. This period is recognized as one of the most crucial for a high final yield of winter wheat [44,45].

Soil Signal
As shown in Figure 2d, the pilot field exhibits a particularly wide range of soil ECa values, between 30 and 130 mS m −1 , which is expected to affect final yield production [46,47]. Indeed, as shown in Figure 10, final yield is highly correlated with soil ECa; higher productivity occurs in areas of lower ECa and vice versa. In an attempt to depict this soil variability through a remote-sensed soil signal, the correlations of Sentinel-2 bands (reflectance) as well as several indices from a single date before the crop establishment (i.e., during bare ground period) were examined against ECa and final yield. As shown in Table 1, among Sentinel-2 bands, band 12 exhibits the highest correlation with ECa, followed by band 11. This result may be ascribed to the fact that both bands incorporate water signals. Accordingly, NMDI-incorporating both 11 and 12 bands-shows the highest correlation with soil conductivity, followed by NDWI which incorporates only band 12. However, it is worth to note that NDVI of bare ground also shows a high correlation with ECa.
Agronomy 2021, 11, x FOR PEER REVIEW 13 of 21 Figure 9. Relationships between water potential (Ψ) and Sentinel-2 water indices NDWI (a) and NMDI (b). Data concern dates with both Sentinel-2 acquisitions and field measurements. Each point corresponds to one Ψ measurement performed in one plot at one date and the mean from the corresponding 9 Sentinel-2 pixels ± standard deviation.
Accordingly, in order to examine the importance of incorporation of a satellite water signal in the model, the correlation of the two water indices against final yield throughout the growing period were investigated. For both indices the most significant water signal (highest seasonal correlation) was found for 29 April (R = 0.641 and 0.550 for NDWI and NMDI, respectively), shortly after flowering. This period is recognized as one of the most crucial for a high final yield of winter wheat [44,45].

Soil Signal
As shown in Figure 2d, the pilot field exhibits a particularly wide range of soil ECa values, between 30 and 130 mS m −1 , which is expected to affect final yield production [46,47]. Indeed, as shown in Figure 10, final yield is highly correlated with soil ECa; higher productivity occurs in areas of lower ECa and vice versa. In an attempt to depict this soil variability through a remote-sensed soil signal, the correlations of Sentinel-2 bands (reflectance) as well as several indices from a single date before the crop establishment (i.e., during bare ground period) were examined against ECa and final yield. As shown in Table 1, among Sentinel-2 bands, band 12 exhibits the highest correlation with ECa, followed by band 11. This result may be ascribed to the fact that both bands incorporate water signals. Accordingly, NMDIincorporating both 11 and 12 bands-shows the highest correlation with soil conductivity, followed by NDWI which incorporates only band 12. However, it is worth to note that NDVI of bare ground also shows a high correlation with ECa.

Model Configuration
The combination of the components described above, i.e., plant, water, and soil signals, was examined against final yield by multilinear regression analysis in the pilot field. Several combinations of VIs corresponding to the three signals were tested and are presented in Table 3. More specifically, concerning plant signal, NDVI and EVI integrals from flowering to the end of the growing period (20 April-31 May) or the best 20 days period after flowering (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25) and their maximum seasonal value (max) were examined. Accordingly, for water signal, NDWI and NMDI for the date with the maximum correlation with final yield (29 April), i.e., the most critical water period, were incorporated. Finally, NDWI and NMDI for a single date before sowing were tested as soil signals.
As shown in Table 3, the best final yield estimation is achieved by a model version combining the NDVI integral from 20 April to 31 May as plant signal, NDWI at 29 April as water signal, and NDWI of bare ground as soil signal (Figure 11a). However, model versions with NMDI as water and/or soil signals show similar high correlations with final yield. For EVI used as plant signal, the best estimation of final yield was also achieved by the 20 April-31 May integral. It is also noteworthy that the performance of a model version with a single date NDVI (Max) for all water and soil signal combinations, was remarkably high further extending the potential for early season yield prediction. Table 3. Coefficients of determination (R 2 ) and root mean square error (RMSE, kg ha −1 ) between measured final yield and various model versions incorporating different plant, water, and soil signals. Best 20 days integral concerns the 6-25 May period. Data concern the pilot field, corresponding to 221 pixels. The most significant correlations are indicated in bold.  Table 3. More specifically, concerning plant signal, NDVI and EVI integrals from flowering to the end of the growing period (20 April-31 May) or the best 20 days period after flowering (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25) and their maximum seasonal value (max) were examined. Accordingly, for water signal, NDWI and NMDI for the date with the maximum correlation with final yield (29 April), i.e., the most critical water period, were incorporated. Finally, NDWI and NMDI for a single date before sowing were tested as soil signals. As shown in Table 3, the best final yield estimation is achieved by a model version combining the NDVI integral from 20 April to 31 May as plant signal, NDWI at 29 April as water signal, and NDWI of bare ground as soil signal (Figure 11a). However, model versions with NMDI as water and/or soil signals show similar high correlations with final yield. For EVI used as plant signal, the best estimation of final yield was also achieved by the 20 April-31 May integral. It is also noteworthy that the performance of a model version with a single date NDVI (Max) for all water and soil signal combinations, was remarkably high further extending the potential for early season yield prediction. Table 3. Coefficients of determination (R 2 ) and root mean square error (RMSE, kg ha −1 ) between measured final yield and various model versions incorporating different plant, water, and soil signals. Best 20 days integral concerns the 6-25 May period. Data concern the pilot field, corresponding to 221 pixels. The most significant correlations are indicated in bold.  Figure 11. Relationships between measured yield and best model estimates based on NDVI (a) or EVI (b) as plant signal. For both models, water and soil signals correspond to NDWI on 29 April and NDWI before sowing, respectively. Data concern the pilot field, corresponding to 221 pixels. The black line corresponds to the 1:1 line and the red is the best-fit line. The equations of the modelled yield obtained from the multiple linear regressions are indicated in the graphs.

Model Extension/Evaluation
In order to extend/evaluate the models described in the previous step for the pilot field, their versions were tested in 31 fields during the 2017-2018 and 2018-2019 growing periods, with a total number of 7436 pixels. It has to be noted that the pilot field was a low productivity field with a yield range between 1000 and 3500 kg ha −1 , whereas the 31 fields used in the extension/evaluation step covered a yield range from 2000 to 7500 kg ha −1 . Accordingly, the full equations determined through multiple linear regression analysis in the pilot field were not used per se in the 31 fields, but only the same signal parameters determined previously were used in a new regression analysis. Therefore, in this step, evaluation concerns the use of the same signal parameters, while extension concerns the new multilinear regressions with theses parameters.
As shown in Table 4, the best-performing versions of the model concern the same plant signals as in the previous step (NDVI and EVI integrals for the 20 April-31 May period), but NMDI instead of NDWI as water and soil signals (R 2 = 0.613, RMSE = 549 for NDVI and R 2 = 0.629, RMSE = 538 for EVI, Figure 12). However, in contrast to the previous step where NDVI was found to perform better than EVI as plant signal, in this step EVI performed better, even though the differences from NDVI was marginal. Accordingly, the performance of the model version with maximum VI as plant signal (max) was very good (R 2 = 0.584, RMSE = 570 for NDVI and R 2 = 0.587, RMSE = 568 for EVI), corroborating its applicability for early and simple yield predictions. Table 4. Coefficients of determination (R 2 ) and root mean square error (RMSE, kg ha −1 ) between measured final yield and various model versions incorporating different plant, water, and soil signals. Best 20 days integral concerns the 13 April-2 May period. Data concern 31 fields during two growing periods, corresponding to 7436 pixels. The most significant correlations are indicated in bold.  For both models, water and soil signals correspond to NDWI on 29 April and NDWI before sowing accordingly. Data concern 31 fields during two growing periods, corresponding to 7436 pixels. The black line corresponds to the 1:1 line and the red is the best-fit line. The equations of the modelled yield obtained from the multiple linear regressions are indicated in the graphs.

Plant, Water, and Soil Signals
In this study, a two-step process for the development of a single model for wheat yield prediction based on Sentinel-2 satellite data was followed. The determination of the For both models, water and soil signals correspond to NDWI on 29 April and NDWI before sowing accordingly. Data concern 31 fields during two growing periods, corresponding to 7436 pixels. The black line corresponds to the 1:1 line and the red is the best-fit line. The equations of the modelled yield obtained from the multiple linear regressions are indicated in the graphs.

Plant, Water, and Soil Signals
In this study, a two-step process for the development of a single model for wheat yield prediction based on Sentinel-2 satellite data was followed. The determination of the satellitederived parameters used in the modelling process was based on field ecophysiological measurements performed in a pilot field. The plant signals obtained from NDVI and EVI showed a high correlation with the field-measured canopy reflectance. Sentinel-2 performance was evaluated on the basis of VIs intercomparisons with the hyperspectral JAZ canopy spectroradiometer of 0.3 nm spectral analysis, which operated very close to the canopy avoiding any atmospheric disturbance effects ( Figure 3). Additionally, both NDVI and EVI well depicted the seasonal development of the winter wheat crop ( Figure 5). Durum wheat in south Europe is planted in late autumn, usually during November or December, seedlings emerge shortly after, and tillering is completed during the winter period. Flowering is considered a critical period of the crop [48,49] and usually occurs in the middle of April in the Thessaly plain. This growth profile was closely followed by both field-measured LAI and satellite VIs, which showed low values during winter, gradually increasing until April flowering period, followed by a rather steep decrease during May seed-ripening period.
In order to identify the optimum period for predicting potential yield, the daily VI data produced by a timeseries interpolation process for each Sentinel-2 pixel were examined against yield maps recorded by the harvester. The results showed that for both NDVI and EVI virtually no correlation was found between the whole growing period integral (1 December-31 May) and yield, while strong correlations appeared when the integral from flowering to the end of the growing period, or the best 20 days integral after the flowering period, was taken into account ( Figure 7). During this particular period, plants start to lose their chlorophylls [50], so VIs were declining ( Figure 5). For the Mediterranean region, this period is considered crucial in determining the final yield because the plants direct their photosynthetic products to the seed [2,49]. The results are in contrast with findings from Ren et al. [22] and Becker-Reshef et al. [9] who declared that the best correlation between NDVI and wheat yield coincides with the period of highest LAI achieved a few days before flowering. However, corroborating our results, Lopresti et al. [24] reported that prediction of wheat yield is best 30 days before harvest, after the stages of heading and flowering. Seggara et al. [51] investigated the optimum period for winter wheat yield estimation in Spain and found that the predictions made on the heading stage outperformed the predictions on tillering or maturing stages, but they did not make any observations during the stage of flowering.
The water signal incorporated in our model was derived from two satellite water indices (NDWI, NMDI), whose efficiency to monitor variations in plant water content was examined by comparing them with field-measured plant water potential. Although NMDI is considered more advantageous compared to NDWI as it incorporates two water bands, it showed a weaker correlation with Ψ (R = 0.775 vs. R = 0.844, respectively, Figure 9). Accordingly, NDWI performed better than NMDI as a water signal in the pilot field, but NMDI outperformed NDWI in the application of the model in the second modelling step involving 31 fields during two growing periods, i.e., covering a wider range of conditions. Soil signal is critical in the construction of a yield prediction model, especially in the case of remarkably variable soil properties, such as ECa. The pilot field of the present study displayed a particularly high range of soil ECa between 30 and 130 mS m −1 , which showed a strong negative linear relationship with wheat yield (Figure 10). Therefore, to address the challenge of capturing such a wide range of ECa through remote-sensed data, the correlations of Sentinel-2 bands and several indices from a single date before the crop establishment (i.e., during bare ground period) were examined against ECa and final yield. NDWI and NMDI of bare ground presented high correlations and were selected for further testing in the procedure of model development. As in the case of the water signal, NDWI performed better than NMDI in the pilot field, and NMDI better in the 31 fields step, indicating its superiority for wide range conditions applications.

Model Development and Validation
Yield prediction of durum wheat has been of major interest both from an economical and managerial point of view. Its accurate prediction requires identification of crop variability even at a within-field scale [8]. The development of remote sensing technologies based mainly on satellite observations made possible the assessment of within-field variability. The Sentinel-2 satellites, operated by the ESA Copernicus program, introduced a new era in open source, high spatial resolution, and frequent earth observations. This advancement is promising not only for capturing within-field variability but also for monitoring of small-sized fields. The latter is crucial for arable lands of the Mediterranean basin due to the small area of farms, the monitoring of which requires high spatial resolution of the remotely sensed data.
The evolution of yield prediction models started with laborious ecophysiological measurements-based simulations and continued with the first generation of satellite-based measurements incorporating low resolution data. During the last decade, complex models have been proposed that require numerous inputs such as soil properties, crop variety information, and meteorological data. Quite often though, the availability or the accessibility of such a broad information is limited. Therefore, there is an essential need for the construction of prediction models that are minimal in inputs and incorporate variables which are easily accessed. In countries like Greece, for example, meteorological data of arable land are not always available or easily accessible. Additionally, between-years climate fluctuation may be unpredictable in the Mediterranean region, especially in the framework of the ongoing climate change. Thus, finding the balance between simplicity and accuracy may result in advantageous next-generation models, with the following characteristics: (i) have as few as possible inputs, without compromising accuracy, (ii) preferentially use satellite data which are easily accessible and present a good correlation with plant function and productivity, (iii) efficiently fit to small farms but could be extended to cover a range of similar environments and cultivation years-high spatial and temporal analysis and extension. In line with this, satellite-based models which utilize the VIs provided by new generation satellites, like ESA's Sentinel-2, are expected to be proved valuable, as they provide information of high spatial and temporal scales.
The modeling procedure presented in this study aimed at addressing the abovementioned needs. The satellite-derived parameters corresponding to plant, water, and soil signals that were determined in the pilot field were further extended/evaluated on 31 other durum wheat fields located in the Thessaly plain throughout two different cultivation years. As shown in Figures 7 and 12, the pilot field was a low-productivity field (yield between 1000 and 3500 kg ha −1 ), whereas the numerous fields used in the extension/evaluation procedure covered a much wider yield range (2000 to 7500 kg ha −1 ). Accordingly, the models developed in the pilot field were not used per se in the other fields, but new regressions with the same signal parameters determined previously were made. Our results show that the overall best performing model version involves the EVI integral for the 20 April-31 May period as a plant signal and NMDI at 29 April and before sowing as water and soil signals, respectively (R 2 = 0.629, RMSE = 538). This model version may be used for yield estimation but not for yield prediction since it incorporates data until the end of the growing period. However, as has been shown (Table 4), the model versions with the maximum seasonal VIs values as plant signal performed almost equally well; the best-performing version involved the max EVI as plant signal and NMDI as water and soil signals (R 2 = 0.587, RMSE = 568). Usually, the maximum seasonal VIs values are recorded during the last ten days of April ( Figure 5), i.e., 30 to 50 days in advance before harvest. Additionally, for this model version only one satellite image is needed for the estimation of the plant signal parameter, significantly simplifying the application of the model. Future work will focus on the optimization of the proposed model through incorporation of additional fields with different soil characteristics and management practices and/or in different areas and cultivation years. Thus, by increasing the conditions range, the equations determined by multilinear regression analysis are expected to become finetuned leading finally to a robust widely applied model. To this end, the determination of the significant signal parameters involved in the modeling approach, which was made in this study, is the first step of the process. Even if it is proved that this approach cannot be applied for widely different conditions, model parameterization for relatively small areas with similar conditions will be a rather simple process. Additionally, the estimation/prediction accuracy of the models presented in this study may be further enhanced if they are used as the basis in productivity models like the ones following the light use efficiency approach, which may better account for climatic and regional variability, increasing however the complexity of the process.

Conclusions
The high spatial and temporal resolution of VIs retrieved from Sentinel-2 imagery enabled the development of several models for durum wheat yield estimation/prediction with good accuracy. EVI was found to function marginally better from NDVI as a plant signal, while NMDI outperformed NDWI as water and soil signals. The best performing model version for yield estimation involved the EVI integral for the 20 April-31 May period as a plant signal and NMDI at 29 April and before sowing as water and soil signals, respectively (R 2 = 0.629, RMSE = 538). Accordingly, the best performing single date model for yield prediction was attained by a version involving the seasonal max EVI as plant signal and NMDI as water and soil signals (R 2 = 0.587, RMSE = 568) 30 to 50 days in advance before harvest. The main advantages of the presented model are its simplicity, the use of easily accessible satellite data of high spatial resolution, and its accuracy in small farms of the Thessaly plain, as was proved in the validation process. Therefore, this durum wheat yield prediction model can be a useful tool for stakeholders, especially dealers, traders, and pasta food companies.