## 4 Discussion  
Testing the applicability of empirical vegetation index models in different years of data allows for evaluating the reliability and reproducibility of that model through time. In this case the stepwise multilinear model split between the training dataset, season 2018, and validation dataset, season 2019, revealed the efficiency of this approach for at least a two-year period in Navarre. This is uncommon, as generally such empirical models are developed and validated with same-year data [16; 21], something that likely happens due to the difficulty to obtain data from various growing seasons, and thus gives robustness to the techniques and analyses that we have presented here. The full operationality of Sentinel-2 constellation reached in 2018 allowed focusing on phenological stages thanks to its improved frequency, as the results presented here show. Until this recent date, phenology-based grain yield prediction was unlikely to be assessed with alike satellites (i.e., Landsat 8). Our results suggest that heading is the most suitable phenological stage in order to study empirical relations of vegetation indices with grain yield for rain-fed wheat in Navarre. This is consistent with other ground spectral measurements for estimating wheat grain yield under rain-fed conditions as reported by [PERSON] et al. [56]. Moreover, this stage was also found optimal for grain yield prediction with proximal remote sensing (i.e., UAV) and equivalent growing conditions [57]. We argue that the results here obtained showed that focusing on phenology for grain yield estimation is now possible using Sentinel-2 imagery and not only applicable with ground and proximal sensing. Furthermore, the combination of climatic (i.e., temperature, growing degree days) and satellite data has resulted in an adequate approach for estimating phenology, as previous studies have suggested [58; 59], and wheat grain yield. The minimum average of pixel values at heading stage of MSAVI and the maximum of RDVI provided the most optimal combination of vegetation indices in the stepwise multilinear regression produced in this study. The improved features of Sentinel-2 regarding spectral and spatial resolution  
Figure 6: GWR local R\({}^{2}\) distribution showed on the four analyzed parameters (altitude, slope, aspect, and rainfall maps, in order from top to bottom) zoomed in to zone II (north). GWR: geographically weighted regression. UTM coordinates are expressed in meters. have allowed increasing the VIs able to be calculated and the most optimal statistic summary, which with coarser spectral and spatial resolutions would be very limited. The improved capacities of these indices could be explained due to the reduced soil background influences, regarding MSAVI, and, as detailed by [PERSON] (2018), the capacity to minimize background reflectance effects of RDVI. In this case, MSAVI together with RDVI show the ability to maintain sensitivity to total vegetation biomass for fully developed canopies corresponding to heading, while other vegetation indices might be saturated at this phenological stage. These vegetation indices have already been successfully correlated with wheat performance using remote sensing data. In the scientific literature, we find examples of RDVI, which is used for green biomass estimation as it is sensitive to vegetation coverage ([PERSON] and [PERSON], 2018) and has also been successfully used in empirical models estimating wheat yields ([PERSON] and [PERSON], 2018). Around the phenological stage of heading, [PERSON] et al. (2018) reported a correlation of 0.79 (R\({}^{2}\)) between RDVI and wheat yield, while [PERSON] et al. (2018) found a correlation of 0.76 (R\({}^{2}\)), which show the relative high correlation ability of RDVI with wheat yield. Examples are also available for MSAVI, which has shown relative high correlation with wheat yields ([PERSON] and [PERSON], 2018; [PERSON] and [PERSON], 2018).