Assessing the Performance of Satellite-Based Models for Crop Yield Estimation in the Canadian Prairies

Abstract

Timely monitoring of crop production using a remote sensing-based approach offers promise toward enhancing food security. Statistical models developed using satellite data typically employ a single vegetation index from a single sensor for yield estimation. With the increasing availability of satellite datasets, there is now an opportunity to investigate the potential of available vegetation indices from different sensors in estimating yields. The key objective of this study was to develop a best-performing satellite-based yield model for the Canadian Prairies for wheat, barley, and canola, trained using municipality-level data from 2009 to 2019. We tested the statistical performance of models built using (a) indices from different sensors (Landsat and Sentinel-2), (b) indices sensitive to different yield properties, and (c) single versus multiple vegetation indices. Results showed Landsat-NDWI as the best performing single-index across all indices and sensors for each crop. Sentinel-2-EVI performed best for wheat and canola and Sentinel-2-SR for barley; but these models were built using only 4 years of data. We found that best-performing single-index models recorded similar predictive accuracy as multi-index models during model validation. The results from this work suggest that satellite-based yield estimation can be improved by selecting the right index related to different crop properties.

RÉSUMÉ

La surveillance en temps voulu de la production agricole à l’aide d’une approche basée sur la télédétection est prometteuse pour l’amélioration de la sécurité alimentaire. Les modèles statistiques développés à partir de données satellitaires utilisent généralement un seul indice de végétation provenant d’un seul capteur pour l’estimation des rendements. Avec la disponibilité croissante des ensembles de données satellitaires, il est désormais possible d’étudier le potentiel des indices de végétation disponibles provenant de différents capteurs pour l’estimation des rendements. L’objectif principal de cette étude était de développer un modèle de rendement basé sur les données satellitaires les plus performantes pour les Prairies canadiennes pour le blé, l’orge et le canola, entraîné à l’aide de données au niveau des municipalités de 2009 à 2019. Nous avons testé la performance statistique des modèles construits en utilisant (a) des indices provenant de différents capteurs (Landsat et Sentinel-2), (b) des indices sensibles à différentes propriétés de rendement, et (c) un seul indice de végétation par rapport à plusieurs. Les résultats ont montré que Landsat-NDWI était l’indice unique le plus performant parmi tous les indices et capteurs pour chaque culture. Sentinel-2-EVI a donné les meilleurs résultats pour le blé et le canola et Sentinel-2-SR pour l’orge; mais ces modèles ont été construits en utilisant seulement 4 années de données. Nous avons constaté que les modèles à indice unique les plus performants ont enregistré une précision prédictive similaire à celle des modèles multi-index lors de la validation du modèle. Les résultats de ce travail suggèrent que l’estimation du rendement par satellite peut être améliorée en sélectionnant le bon indice lié aux différentes propriétés des cultures.

Introduction

Timely, consistent, and reliable estimation of crop yields are key to making strategic decisions related to food security planning from local to global scales, especially under changing growing conditions. Crop yield is influenced by interactions between multiple factors such as weather, soil health, and management choices (e.g., irrigation, fertilizers). The influence of different factors that drive yields are reflected in crop growth conditions that can be directly measured using remotely sensed satellite imagery. Satellite data provides near-real-time intelligence on vegetation health status suitable for continuously monitoring and assessing crop yields over space and time (Sishodia et al. 2020). Historical and near-real time crop production insights are critical inputs for decision making by different agricultural stakeholders (farmers, crop insurers, policy makers, commodity traders) (Weiss et al. 2020; Karthikeyan et al. 2020).

Different agricultural stakeholders (e.g., producers, commodity traders, government, crop insurers) require different pieces of intelligence on agricultural production toward achieving diverse targets (Table 1). For example, producers want to optimize production quantity using historical and near-real time data, crop insurers aim to design satellite-based insurance contracts, and government or policy makers are interested in strengthening food security using satellite-based food monitoring systems. Accordingly, the role of satellite vegetation indices in providing context-specific yield insights in addition to accurate yield estimates is critical. Currently, most agricultural stakeholders rely on models using traditional single indices such as NDVI and EVI which may not be well suited to monitoring nutrient-stressed crops (e.g., in nutrient-limited regions) or water-stressed crops (e.g., in irrigation-limited regions). The careful selection and inclusion of additional information by combining multiple, different indices measuring different crop properties could be useful for supporting existing satellite-based agricultural intelligence capabilities.

Table 1. A selection of use-case applications of satellite-based vegetation indices by different agricultural stakeholders.

Stakeholder	Use case	Problem description	Spatial scale	Temporal scale	Sensor	Indices	Reference set
Farmer	Precision farming	Optimize yields	Field/Farm	NRT	Sentinel-2	NDVI, Red-edge VIs	A
		Track water stress	Field/Farm	NRT	Sentinel-2, Landsat	NDWI, NDVI, EVI	B
		Monitor soil health	Field/Farm	Both	Landsat	NDVI, NDWI	C
		Track nutrient stress	Field/Farm	NRT	Skysat, Landsat	GCVI, GNDVI	D
		Monitor disease	Field/Farm	NRT	Sentinel-2	Red-edge VIs	E
	Yield forecasting	Optimize yields	Field/Farm	NRT	MODIS, Sentinel-2	EVI2, NDVI	F
	Yield mapping	Build historical yield database	Field/Farm	HDA	Sentinel-2, Landsat	NDVI, GCVI	G
Crop insurers	Risk management	Design index-based insurance	Sub-region, region	HDA	AVHRR, MODIS, Landsat	NDVI, NDWI	H
	Insurance planning	Regional yield disaggregation to local scale	Sub-region	HDA	MODIS	EVI	I
Government	Strengthen food security	Forecast crop yields	Region, national	NRT	MODIS, Sentinel-2	NDVI, EVI, Red-edge VIs	J
		Map croplands	Region, national	HDA	MODIS, AVHRR	NDVI, NDWI	K
Commodity traders	Yield forecasting	Commodity arrival prediction	National, global	Both	MODIS	NDVI	L
	Assess food sustainability scenarios	Map commodity lands	National	HDA	MODIS, AVHRR	NDVI, EVI	M
International Org.	Monitor food systems	Forecast crop yields	National, global	NRT	MODIS	NDVI, NDWI	N
		Develop food database	National, global	HDA	Sentinel-2	NDVI	O
		Assess ecosystem services	National, global	HDA	Landsat	NDVI	P

NRT is near-real time and HDA is historical data archive.

Reference studies are listed below.

A: Nguyen et al. 2022; Hunt et al. 2019; Zhao et al. 2020; Cavalaris et al. 2021.

B: Potgieter et al. 2021; Pantaleoni et al. 2007; Cavalaris et al. 2021; Cicek et al. 2010; Soriano-González et al. 2022.

C: Satir and Berberoglu 2016; Tadesse et al. 2021.

D: Burke and Lobell 2017; Franz et al. 2020; Kayad et al. 2019.

E: L. Liu et al. 2020; Zheng et al. 2018; Yang 2020.

F: Bolton and Friedl 2013; Nazir et al. 2021.

G: Hunt et al. 2019; Franz et al. 2020.

H: Kölle et al. 2022; Möllmann et al. 2020; Bokusheva et al. 2016; Masiza et al. 2022; Borgogno-Mondino et al. 2019; Makaudze and Miranda 2010; Turvey and Mclaurin 2012; Benami et al. 2021; Ha et al. 2022.

I: Shirsath et al. 2020.

J: Chipanshi et al. 2015; Qader et al. 2018; Karlson et al. 2020.

K: Becker-Reshef et al. 2010; Forkuor et al. 2017; Salmon et al. 2015; Beneke et al. 2019.

L: Prasad et al. 2019; Doraiswamy et al., n.d.

M: Condro et al. 2020; Brown and Shahriar Pervez 2014.

N: Becker-Reshef et al. 2010.

O: Bontemps et al. 2015.

P: Shi et al. 2020.

Most studies using a satellite-based approach for estimating crop yields employ a single vegetation index for deriving yield estimates (Liu et al. 2020; Bolton and Friedl 2013; Wall et al. 2008; Burke and Lobell 2017; Kouadio et al. 2014). Traditional vegetation indices derived from popular satellite data sources (AVHRR¹, MODIS², Landsat) include normalized difference vegetation index (NDVI) and two-band enhanced vegetation index (EVI2), both tracking the biomass component of yield. MODIS-based EVI-2 has been found to perform better than NDVI at crop yield prediction previously (Liu et al. 2020). Some studies have found normalized water difference index (NDWI), which is sensitive to crop water content, to be an effective predictor of yields in more arid regions (Bolton and Friedl 2013). Other studies note that a green chlorophyll vegetation index (GCVI), tracking variation in chlorophyll concentration is better suited than NDVI or EVI to capture the effects of nutrient stress on resulting yields (Burke and Lobell 2017; Azzari et al. 2017). GCVI is used because it has been suggested to approximate LAI well (Nguy-Robertson et al. 2012). Another study using combined Landsat and Sentinel-2 data reports that best results were obtained by a model with NIR and red spectral bands (Skakun et al. 2019). Marszalek et al. (2022) highlights the importance of indices tracking crop water and nitrogen content derived from red-edge and near-infrared bands of the Sentinel-2 sensor for yield prediction at the field level.

Few previous studies have evaluated the potential of combining multiple trait-specific vegetation indices for estimating yields or compared different satellite data sources and their value in estimating yields for multiple crops. Most previous studies (with a few exceptions; Zhao et al. (2020), Ulfa et al. (2022), Cavalaris et al. (2021) have focused on studying one crop type (region-specific major crops) using a single satellite data source. Zhao et al. (2020) investigated the potential of combining biomass and chlorophyll indices for wheat yield prediction at the field level using Sentinel-2 data. They found a small increase in model predictive skill when different indices at their respective peaks were combined. They recommend further testing of a multi-index modeling approach in different regions with diverse crop management practices. Similar results were also noted by Ulfa et al. (2022) in modeling field-level wheat yields using multi-indices derived from Landsat data. Currently, a multi-index crop yield estimation approach remains mostly constrained to the field level and one crop type with limited data records.

In this study, we investigate the statistical performance of satellite-based crop yield models including single vegetation indices and multiple vegetation indices tracking different yield traits derived using different satellite datasets. We train and cross-validate all models using observed historical yield data (2009–2019) at the municipality level for three key crops (wheat, barley, canola) in the Canadian Prairies. The main objective of this study was to develop a best performing satellite-based crop yield model. Our study had the following three sub-objectives: (1) compare how well different vegetation indices available from different satellite sensors predict yields, (2) compare accuracy of single-index versus multi-index models for each sensor, and (3) evaluate if model performance improved when different indices sensitive to different yield properties were included in the same model.

Materials and methods

Study area: Canadian Prairies

The Canadian Prairies are an economically important region of Canada with field crop production accounting for almost half of agricultural operations in the region. Spring wheat, barley, and canola are the largest commodities grown by planted acreage. The study area covers the three provinces of Alberta, Saskatchewan, and Manitoba (Figure 1). Regional agro-climatic conditions vary year-to-year resulting in differing growing environments. Precipitation is highly variable (300–500 mm per annum), both spatially and temporally, leading to moisture deficit and crop health stress in irrigation-limited areas. The growing season for annual spring crops begins with seeding in April and concludes with harvesting between early to late fall. Majority of the agricultural land in Manitoba is located within the moist black and dark gray soil zones. Agricultural regions of Alberta and Saskatchewan are located across differing soil zones experiencing varying interactions between soil and climate resulting in diverse crop growth environments.

Figure 1. Spatial distribution of sensor-specific selected mean vegetation indices observed at the peak of the growing season at municipality level of the Canadian Prairies for wheat, barley and canola during the 2016–2019 period. Regions in white represent the absence of cropland within a given municipality.

Crop yield and satellite data (MODIS, Landsat, Sentinel-2)

There are a total of 503 municipalities within the Canadian Prairie Provinces. Municipality-level yield data for spring wheat, barley, and canola from 2009 to 2019 were obtained from different sources. Crop yield datasets for Saskatchewan³ and Manitoba⁴ were accessed from open-source databases maintained by provincial governments. Yield data for Alberta was provided by the Agriculture Financial Services Corporation (AFSC). Crop yield in this study refers to the end-of-season yield recorded by each reporting agency. As crop-specific harvested area at the municipality level is publicly unavailable, we were unable to set an area-related inclusion or exclusion threshold for municipalities. However, as we are modeling yields, estimated as production per unit harvested area, the outcome variable is normalized by area.

The sources of satellite imagery used in this study are MODIS (Terra) at 250 m spatial resolution, Landsat (7&8) at 30 m spatial resolution, and Sentinel-2 (A&B) at 10 m spatial resolution. MODIS and Landsat data records are available from 2009 to 2019. As Sentinel-2 data becomes available from 2016 onwards, Sentinel-2 data covers the period 2016–19. All satellite data products used were level-2 orthorectified, surface reflectance products. The satellite datasets provide contrast in temporal and spatial resolution that could affect the performance of satellite-based yield estimation methods. MODIS and Sentinel-2 imagery have a finer temporal resolution with daily and near-daily (every 3 d) revisit cycles respectively, compared to Landsat (7&8) imagery which is available every 16 d. As a previous study in the Prairies found seasonal peak EVI2 to perform better at crop yield estimation than other metrics derived using MODIS data (Liu et al. 2020), we extracted a single-index, EVI2, using MODIS data. With Landsat and Sentinel-2 data, we retrieved seven and ten trait-specific vegetation indices respectively (Table 2). Annual crop-specific masks freely available from the Canadian annual crop inventory at 30 m spatial resolution from 2011 onward (crop masks are available at 56 m spatial resolution for years 2009 and 2010) were applied prior to calculation of crop-specific indices. The crop-specific masks captured both spatial and temporal cropping conditions of only cropping areas within each municipality targeting the specificity of data used for model training toward improving the quality of estimates. A sensor-specific cloud mask with a threshold of < 10% cloud coverage was applied to all satellite imagery to obtain cloud-free pixels prior to calculation of indices. All suitable Landsat images were selected to derive maximum monthly composites toward the calculation of indices. Due to a 16-day repeat cycle associated with the Landsat sensor, there were missing values for some municipalities. MODIS- and Sentinel-2-based indices were extracted as maximum weekly composites. A previous study noted that yield estimation using Sentinel-2 data extracted at a weekly resolution had a superior performance than monthly resolution data as a weekly temporal scale better reflected the crop growth status (Marszalek et al. 2022). With each image collection, vegetation indices were extracted as a spatial mean at the municipality level. All satellite datasets were accessed and processed in Google Earth Engine (GEE) due to the platform’s capability in handling large volumes of satellite imagery (Gorelick et al. 2017). A selected summary of spatial distribution of satellite vegetation indices at municipality level used as model inputs is shown in Figure 1.

Table 2. Sensor-specific vegetation indices tracking different properties of crop yield.

Trait tracked	Index	Sensor	Ref. study
Biomass	EVI2	MODIS	(Liu et al. 2020)
	NDVI, EVI	Landsat, Sentinel-2	(Zhao et al. 2020; Ulfa et al. 2022)
Water	NDWI, GI	Landsat, Sentinel-2	(Bolton and Friedl 2013; Deines et al. 2017)
Chlorophyll	GCVI, GNDVI, SR	Landsat, Sentinel-2	(Burke and Lobell 2017; Kayad et al. 2019)
	NDRE1, NDRE2, TCARI	Sentinel-2	(Zhao et al. 2020)

Crop yield estimation

We use a statistical linear regression modeling approach common in the literature to estimate yields from satellite indices (Johnson et al. 2021; Basso and Liu 2019). Two yield models, single-index and multi-index, for each crop were developed. Crop-specific biomass index, water index and chlorophyll index were combined in a multi-index modeling framework (Equation 1(1) $$y=\mathrm{ }{\beta }_0+{\beta }_{1}BI+\mathrm{ }{\beta }_2\mathrm{ }WI+{\beta }_{3}CI+\mathrm{ }{\beta }_{4}M+{\beta }_{5}T$$(1) ) to test for additional improvements over a single-index model (Equation 2(2) $$y=\mathrm{ }{\beta }_0+{\beta }_{1}I+\mathrm{ }{\beta }_{2}M+{\beta }_{3}T$$(2) ). In each multi-index model, we only included a single-index representing a specific yield trait as a potential predictor to avoid multi-collinearity issues. We also conducted additional testing for multicollinearity between different indices using Variance Inflation Factor (VIF) before including respective indices in multi-index models and found that combinations of indices that had lower VIF values also had lower Akaike information criterion (AIC) values. As we wanted to conduct a comprehensive analysis using a multi-index approach, we tested all possible combinations of trait-specific indices with Landsat data (12 trait-specific index combinations) and Sentinel-2 data (24 trait-specific index combinations) to identify the best combination of indices for yield estimation (Table A7). Models using a single index (Equation 2(2) $$y=\mathrm{ }{\beta }_0+{\beta }_{1}I+\mathrm{ }{\beta }_{2}M+{\beta }_{3}T$$(2) ) for all Landsat- and Sentinel-2-derived indices individually were compared to multi-index models. (1) $$y=\mathrm{ }{\beta }_0+{\beta }_{1}BI+\mathrm{ }{\beta }_2\mathrm{ }WI+{\beta }_{3}CI+\mathrm{ }{\beta }_{4}M+{\beta }_{5}T$$(1) (2) $$y=\mathrm{ }{\beta }_0+{\beta }_{1}I+\mathrm{ }{\beta }_{2}M+{\beta }_{3}T$$(2)

Where I = {B, W, C}

In the regression Equations (1)(1) $$y=\mathrm{ }{\beta }_0+{\beta }_{1}BI+\mathrm{ }{\beta }_2\mathrm{ }WI+{\beta }_{3}CI+\mathrm{ }{\beta }_{4}M+{\beta }_{5}T$$(1) and (2)(2) $$y=\mathrm{ }{\beta }_0+{\beta }_{1}I+\mathrm{ }{\beta }_{2}M+{\beta }_{3}T$$(2) , $y$ represents crop yield, I represents satellite vegetation index, M represents random municipality effects and T represents a linear year trend. BI, WI and CI represent a biomass, water and chlorophyll index respectively.

As imagery from different satellite sources are available at different times during the year, we developed both multi-index and single-index models at different time steps specific to the sensor to gauge which model performed best. With Landsat data, we developed monthly models within each growing season using monthly peak indices. Weekly models were developed for each of twenty-two weeks within each growing season using Sentinel-2 weekly index composites. With MODIS data, we developed a single model using seasonal peak vegetation index as peak stage was found to be the optimal time to obtain yield estimates by a previous study in the Prairies (Liu et al. 2020). We also investigated the performance of multi-stage⁵ against single stage models and found that a single stage model performed better than a multi-stage model for each crop at yield estimation (Figure A6). All models include random municipality effects⁶ to account for local variations at the municipality level due to differing management practices (e.g., irrigation, tillage, cultivars etc.) and linear year trend to account for the effects of technological developments on crop yield.

Model selection and assessments

We used the Akaike information criterion (AIC) as the model selection metric to identify the best-performing multi-index model combination over all possible index combinations considered. AIC as a model selection metric weight both model performance and model complexity. For models with the lowest AIC (AICc for small samples) value (Tables A4–A6), we performed leave-one-year-out cross-validation (LOOCV) to estimate the prediction error of chosen models at the municipality level. LOOCV as an out-of-sample validation test allows to evaluate the generalizability of a model in different years with different climatic conditions (Liu et al. 2020; Kouadio et al. 2014). Finally, we compared the prediction accuracy of chosen models using the following metrics: root-mean-square-error (RMSE) measured in tons per hectare, relative root-mean-square-error (RRMSE, %) to mean observed yield, and mean absolute error (MAE) measured in tons per hectare. We used an additional metric, Wasserstein distances (WD), to further evaluate model performance. WD is a distance metric that quantifies the difference or distance between two data distributions (Muskulus and Verduyn-Lunel 2011). The value of WD is equal to 0 if distributions are identical with the value of WD becoming larger as distributions become different. WD was calculated between the data distributions of observed and estimated yields.

Results

Overall performance of multi-index and single-index models at the municipality level

The best performing single-index and multi-index models using Landsat data during 2009–2019 performed quite similarly on cross-validation metrics (RRMSE range: 14.6–19.7% and 14.3–19.4% for single-index and multi-index models respectively across crops) (Table 3). We observed low WD values for best performing crop-specific single-index and multi-index Landsat models showing a good agreement between observed yields and estimated yields. For wheat, the multi-index combination EVI + SR + NDWI showed only a marginal increase in predictive power over NDWI, the best single-index predictor for wheat. A single-index model using NDWI performed slightly better than a multi-index model for barley and canola. For barley and canola, the prediction error of single-index NDWI model (multi-index model) is 18.69% (19.42%) and 16.45% (16.49%) respectively. The best performing multi-index combination for barley and canola were found to be EVI + SR + GI and EVI + SR + NDWI respectively. Overall, we note that a single index tracking water content in crops, NDWI, performed better than all other trait-specific Landsat indices tested for each crop. For example, compared to traditional indices such as NDVI derived using Landsat data, NDWI showed a slightly better performance. The prediction error (RRMSE (%)) of NDWI (NDVI) single-index model is 14.57% (15.14%) for wheat, 18.69% (19.69%) for barley and 16.45% (17.65%) for canola. For modeling cases tested using Landsat data, satellite predictors derived at the height of July (mid-season coinciding with peak of crop growth stage) was observed to perform better at yield estimation than predictors derived during other periods (early- and late-season) within the growing season. Further we note that yield estimation using Landsat data demonstrated improvements over yield estimates derived using MODIS data for each crop. The prediction error (RRMSE (%)) of the best performing Landsat (MODIS) model is 14.30% (16.09%) for wheat, 18.69% (20.16%) for barley and 16.45% (17.47%) for canola. The WD values of crop-specific best performing Landsat model are also low and comparable to WD values of a MODIS model highlighting a good agreement between observed and estimated yields across space and time. Table 3 lists all cross-validation modeling results for the period 2009–2019.

Table 3. Cross validation performance metrics of best performing single-index and multi-index models for full time series (2009–2019). The best performing model (in bold) across model types for each crop was chosen using RMSE (RRMSE) and corresponding $R^2$ and WD values are bolded for the same model.

	Wheat					Barley					Canola
	M	L	L	L	L	M	L	L	L	L	M	L	L	L	L
	EVI2	NDVI	EVI	SI	MI	EVI2	NDVI	EVI	SI	MI	EVI2	NDVI	EVI	SI	MI
RMSE	0.479	0.451	0.437	0.434	0.426	0.814	0.795	0.777	0.755	0.784	0.416	0.419	0.415	0.392	0.392
RRMSE	16.09	15.14	14.67	14.57	14.30	20.16	19.69	19.23	18.69	19.42	17.47	17.65	17.44	16.45	16.49
MAE	0.374	0.359	0.345	0.344	0.337	0.644	0.638	0.626	0.602	0.627	0.337	0.345	0.339	0.317	0.320
$R^2$	0.66	0.64	0.62	0.64	0.61	0.49	0.40	0.40	0.42	0.36	0.49	0.35	0.34	0.40	0.39
WD	0.133	0.142	0.130	0.124	0.127	0.292	0.329	0.309	0.294	0.310	0.147	0.171	0.166	0.146	0.141

M: MODIS-based model and L: Landsat-based model.

SI: best performing single-index model and MI: best performing multi-index model.

SI (L) model: NDWI for wheat, barley and canola.

MI (L) model: EVI + SR + NDWI for wheat, EVI + SR + GI for barley, EVI + SR + NDWI for canola.

Root-mean-square-error (RMSE) in t/ha, RMSE relative to mean observed yield (RRMSE (%)), mean absolute error (MAE) in t/ha, coefficient of determination (R²), and Wasserstein Distance (WD).

We further tested model performances during 2016–2019, matching Sentinel-2 data availability. The length of MODIS and Landsat data records used for model training was shortened to match Sentinel-2 records. During this shorter modeling period, Landsat-based single index GI (RRMSE: 12.79%) and NDWI (RRMSE: 15.02%) for wheat and barley respectively demonstrated good performance over MODIS-based and Sentinel-2-based models (Tables A1 and A2). For wheat and barley, Landsat predictors derived at the height of July remained the optimal time period for obtaining yield estimates. Seasonal peak MODIS-EVI2 (RRMSE: 13.52%) had the lowest prediction error for canola (Table A3). With Sentinel-2 data, we observed that multi-index models were unable to demonstrate huge improvements over single-index models for each crop (Table 4). A single-index model using EVI (RRMSE: 18.17%) tracking biomass content had the lowest prediction error compared to a multi-index (EVI + TCARI + GI) model (RRMSE: 21.35%) for wheat. For barley and canola, a multi-index model performed marginally better than a single-index model. The best performing multi- index combination for barley and canola were found to be EVI + SR + NDWI and EVI + TCARI + GI respectively. The prediction error (RRMSE (%)) of the best performing multi-index (single-index) model is 18.57% (18.88%) for barley and 15.19% (15.31%) for canola. For models using Sentinel-2 data, indices derived toward end of June (week 25) for wheat and barley and at the beginning of August (week 32) for canola had lowest prediction errors. The cross-validation modeling results for the best performing single-index and multi-index models using Landsat and Sentinel-2 data during the 2016–19 period are presented in Table 4. Detailed results for 2016–19 period for wheat, barley and canola are presented in Tables A1–A3 respectively in the appendix. Figures 2 and 3 show the relationship between observed yield and yield estimated by the best performing single-index and multi-index models using Landsat and Sentinel-2 data and by a model using MODIS-EVI2 for periods 2009–2019 and 2016–2019 respectively. We observed an overall pattern of underestimation of high yields and overestimation of low yields using different satellite-based yield models. The spread of data points associated with the MODIS model for each crop is higher compared to Landsat and Sentinel-2 models during both modeling periods.

Figure 2. Benchmarking of best performing single-index and multi-index models using Landsat (7&8) data against MODIS single-index model for full time series (2009–2019) for wheat (Top panel), barley (Middle panel) and canola (Bottom panel). Observed yield shown along the y-axis range from 0 to 8 t/ha for wheat (A–C) and barley (D–F) and from 0 to 4 t/ha for canola (G–I).

Figure 3. Benchmarking of best performing single-index and multi-index models using Landsat (7&8) and Sentinel-2 data against MODIS single-index model for the time period 2016–2019 for wheat (Top panel), barley (Middle panel) and canola (Bottom panel). Observed yield shown along the y-axis range from 0 to 6 t/ha for wheat (A–E), from 0 to 8 t/ha for barley (F–J) and from 0 to 4 t/ha for canola (K–O).

Table 4. Cross validation performance metrics of best performing single-index and multi-index during the time period corresponding to Sentinel-2 data availability (2016–2019).

	Wheat				Barley				Canola
	L	L	S	S	L	L	S	S	L	L	S	S
	SI	MI	SI	MI	SI	MI	SI	MI	SI	MI	SI	MI
RMSE	0.426	0.434	0.605	0.711	0.677	0.718	0.851	0.837	0.401	0.416	0.406	0.402
RRMSE	12.79	13.05	18.17	21.35	15.02	15.91	18.88	18.57	15.12	15.69	15.31	15.19
MAE	0.334	0.339	0.482	0.577	0.529	0.551	0.688	0.685	0.313	0.326	0.328	0.317
$R^2$	0.59	0.57	0.39	0.34	0.49	0.42	0.28	0.29	0.41	0.37	0.37	0.39
WD	0.116	0.120	0.228	0.314	0.215	0.258	0.439	0.456	0.154	0.162	0.133	0.127

L: Landsat-based model; S: Sentinel-2-based model; SI: best performing single-index model; MI: best performing multi-index model.

SI (L) model: GI for wheat, NDWI for barley and canola.

SI (S) model: EVI for wheat and canola, SR for barley.

MI (L) model: EVI + SR + GI for wheat, EVI + SR + NDWI for barley and canola.

MI (S) model: EVI + TCARI + GI for wheat and canola, EVI + SR + NDWI for barley.

Root-mean-square-error (RMSE) in t/ha, RMSE relative to mean observed yield (RRMSE (%)), mean absolute error (MAE) in t/ha, coefficient of determination (R²), and Wasserstein Distance (WD).

Interannual variation in estimated yields

The best performing multi-index and single-index models using Landsat and Sentinel-2 data and MODIS single-index model were used to estimate crop yields at the municipality level over the study period. Figure 4 shows the mean interannual variation in model-estimated yield against observed yield for each crop. Crop yields have increased over the 11-year period with some interannual variation. Yield models based on MODIS and Landsat data were able to capture the overall observed trend over 2009–2019. We observed overestimation of low yields and underestimation of high yields by all models. Landsat-based multi-index and single-index models for canola demonstrated nearly identical estimation performance during the entire study period (Figure 4H). We noted some minor differences in estimation performance of both Landsat-based models for years 2013–2019 and 2009, 2015–2019 for wheat and barley respectively (Figures 4B and 4E). Interannual yield estimates varied by sensor-specific index predictors used to train models over the full period (Figure A1). Sentinel-2 multi-index and single-index model performance during the 2016–2019 period varied by crop type (Figure A2). A marginal difference in performance between both model types was seen for wheat and barley during this short time period (Figures 4C and 4F) whereas for canola, both models exhibited similar performance (Figure 4I).

Figure 4. Mean annual variation in estimated yield against observed yield of wheat (Top panel), barley (Middle panel) and canola (Bottom panel) using MODIS baseline model (A,D,G) and best performing Landsat (B,E,H) and Sentinel-2 (C,F,I) models. Crop yield shown along the y-axis range from 2 to 4 t/ha for wheat (A–C), from 3 to 5 t/ha for barley (D–F) and from 1 to 3 t/ha for canola (G–I). A simple yield average across all municipalities was calculated to get an annual mean yield value for inter-year comparison of model estimated yield to observed yield.

Evaluation of model performance at municipality level

The comparison of model performance over the 2009–2019 period showed that results were quite similar when using a dataset including either the best set of Landsat-based indices or a MODIS index. (Figure 5). The ranges of RRMSE for wheat, barley and canola using Landsat (MODIS) data were −67% (−80%) to 56% (58%), −71% (−93%) to 87% (92%) and −61% (−68%) to 85% (85%) respectively. The overall spatial distribution of yields estimated using different satellite index predictors during periods, 2009 and 2019, for wheat, barley and canola is shown in Figures 6–8 respectively. The spatial distribution of estimated yields for each crop showed that model performances varied slightly across the agricultural municipalities by the index predictor(s) used for model training. The spatial distribution of model estimation error (underestimation, overestimation, no change) for wheat, barley and canola is presented in Figures A7–A9 respectively in the appendix. We note that models were able to capture the increase in yields observed between 2009 and 2019 in some municipalities. Overall, satellite-based yield models were able to satisfactorily reproduce spatial yield variations observed at the municipality level as documented by different yield reporting agencies which is highlighted by low WD values (Tables 3 and 4 and Tables A1–A3). WD values of best performing single-index and multi-index models using Landsat and Sentinel-2 data range between 0.11 − 0.33 and 0.10 − 0.46 respectively.

Figure 5. Ranges of cross validation model prediction error (RMSE) relative to mean observed yield at the municipality level, 2009–2019 period, for best performing model using Landsat data (EVI + SR + NDWI for wheat, NDWI for barley and canola) and MODIS data (EVI2) for each crop.

Figure 6. Spatial distribution of observed yield and estimated yield of wheat using different satellite index predictors for 2009 and 2019. Landsat-NDWI is the best performing single-index. Best performing multi-index combination using Landsat data for wheat is EVI + SR + NDWI. Regions in white represent absence of cropland within a given municipality. Regions in grey represent municipalities for which an model-specific yield estimate is unavailable.

Figure 7. Spatial distribution of observed yield and estimated yield of barley using different satellite index predictors for 2009 and 2019. Landsat-NDWI is the best performing single-index. Best performing multi-index combination using Landsat data for wheat is EVI + SR + GI. Regions in white represent absence of cropland within a given municipality. Regions in grey represent municipalities for which an model-specific yield estimate is unavailable.

Figure 8. Spatial distribution of observed yield and estimated yield of canola using different satellite index predictors for 2009 and 2019. Landsat-NDWI is the best performing single-index. Best performing multi-index combination using Landsat data for wheat is EVI + SR + NDWI. Regions in white represent absence of cropland within a given municipality. Regions in grey represent municipalities for which an model-specific yield estimate is unavailable.

Discussion

This study investigated the potential of combining multiple, different trait-specific satellite-based vegetation indices over using single indices for crop yield estimation using different satellite datasets. Modeling results show that multi-index models based on Landsat and Sentinel-2 data were unable to demonstrate large improvements over models using sensor-specific single indices at yield prediction (Tables 3 and 4). Combining vegetation indices tracking different yield traits namely biomass, water content and chlorophyll content does not seem to provide additional information on crop health and it is likely that overlapping signals results in a lack of improvement in multi-index models. For example, Figure 9 shows the relationship between best performing Landsat-based indices identified for multi-index model (EVI + SR + NDWI) for wheat. High correlations (> 0.75) are observed between biomass index (EVI) and chlorophyll index (SR) (Figure 9A) and between biomass index (EVI) and water index (NDWI) (Figure 9B) (Detailed correlation results for wheat, barley, and canola are shown in Figures A3–A5 respectively in the appendix). Additional research should consider the inclusion of raw band values of satellite images that are used to calculate indices as ratios of bands to investigate if yield estimation quality improves for a multi-index model trained with specific satellite bands. For example, one of the bands included in all chosen indices in the Landsat-based multi-index model for wheat, EVI + SR + NDWI, covers the near-infrared (NIR) spectrum. The NIR-based band tracks the chlorophyll aspect of yield quality and the importance of this band was previously reported (Skakun et al. 2019; Marszalek et al. 2022; Prey and Schmidhalter 2019). Marszalek et al. (2022) notes that Sentinel-2-based raw bands covering the red and NIR spectrum achieved higher accuracies than the index-based approach at yield estimation.

Figure 9. Correlation between best performing Landsat-based vegetation indices used in a multi-index model for wheat (EVI + SR + NDWI).

In this study, single-index models showed prediction performance equivalent to the performance of multi-index models. Model tests during the full study period (2009–2019) showed that Landsat-based single-index models using NDWI performed slightly better than models using multi-indices for barley and canola. For wheat, the prediction performance of the multi-index model was marginally better than the single-index model. A summary of the best-performing Landsat-based single-index and multi-index models during the 2009–2019 period is presented in Tables 5 and 6 respectively. Previous studies have reported similar findings when combining multiple different indices sensitive to different yield traits for prediction of crop yields at the field level (Zhao et al. 2020; Ulfa et al. 2022; Cavalaris et al. 2021; Li et al. 2022). Ulfa et al. 2022 assessed different combinations of Landsat-based biomass and chlorophyll indices and found that a combination of EVI + GCVI performed marginally better than single indices for predicting field level wheat yields. Using Sentinel-2 derived indices, statistical models driven by index combinations such as NDVI + NDRE1, NDVI + NDWI, EVI + NDWI and GCVI + EVI2 have produced slightly better yield estimates over models using a single index (Zhao et al. 2020; Cavalaris et al. 2021; Li et al. 2022). Landsat-based single-index models again demonstrated better prediction performance than multi-index models for each crop during the 2016–19 period. Using Sentinel-2 indices, single-index model using EVI outperformed the multi-index model for wheat whereas prediction performance of single-index models and multi-index models was comparable for barley and canola. As EVI is less prone to saturation with dense crop biomass and better reflects crop health, Sentinel-2-EVI was similarly found to be an important predictor of crop yield in a previous study (Pejak et al. 2022; Li et al. 2022).

Table 5. Best performing crop-specific, sensor-specific single indices.

	2009–2019	2016–2019	2016–2019
	Landsat (7&8)	Landsat (7&8)	Sentinel-2
Wheat	NDWI	GI	EVI
Barley	NDWI	NDWI	SR
Canola	NDWI	NDWI	EVI

Modeling results show that a single-index modeling approach has the potential to perform well compared to a multi-index approach as long as a user employs an appropriate vegetation index suited to the particular satellite product in a yield estimation model. With a single-index approach, identifying the best crop-specific, region-specific and sensor-specific vegetation index is critical for obtaining high quality yield estimates. In contrast, a multi-index modeling approach is beneficial when a user wants to phase out the index identification and selection process. Using a multi-index approach, an user can let the model automatically select and use the index that best tracks a particular yield trait through a growing season (e.g., drought vs non-drought year) to produce estimates. Further testing of a multi-index modeling approach could provide higher utility at yield estimation in operational applications (e.g., forecasting yields, commodity arrival prediction) monitoring diverse agricultural regions. Users deploying a Landsat-based model in an operational environment might observe missing yield estimates in some agricultural regions due to less frequent (16-day) repeat cycle of the sensor and cloud cover. In the future, the application of a Sentinel-2-based model in such a scenario is promising as it has a better repeat cycle (near-daily, that is, every three days) than Landsat. A summary of the best performing crop- specific, satellite sensor-specific single-index and multi-index combinations identified during the two periods tested in this study is presented in Tables 5 and 6 respectively.

Table 6. Best performing crop-specific, sensor-specific multi-index combinations.

	2009–2019	2016–2019	2016–2019
	Landsat (7&8)	Landsat (7&8)	Sentinel-2
Wheat	EVI + SR + NDWI	EVI + SR + GI	EVI + TCARI + GI
Barley	EVI + SR + GI	EVI + SR + NDWI	EVI + SR + NDWI
Canola	EVI + SR + NDWI	EVI + SR + NDWI	EVI + TCARI + GI

Results further revealed the importance of some vegetation indices that could be useful toward improving yield estimates in the Prairie agricultural region. Landsat-based single-index, NDWI, which is sensitive to crop moisture content demonstrated predictive skill similar to multi-index models when tested during the full study period. This finding highlights the capability of NDWI in tracking plant water fluctuations and its resulting impacts on final crop yield. NDWI has been similarly found to be a more effective predictor of crop yields compared to indices such as NDVI and EVI in arid agricultural regions and in croplands with highly variable soil properties (e.g., salinity) (Bolton and Friedl 2013; Cavalaris et al. 2021). This suggests that the quality of existing regional yield forecast models could be improved by using NDWI instead of NDVI or EVI as commonly practiced. For example, the Canadian crop yield forecaster tool currently uses NDVI with climate data as inputs to forecast crop yields in agricultural regions of Canada (Chipanshi et al. 2015). Further exploration of predictors such as NDWI that better represent crop health status over the growing season may be useful for analyzing the spatial variation of crop yields specially in those regions where NDVI does not perform well. With Sentinel-2 derived indices, single-index EVI achieved prediction skill which was close to a multi-index model for wheat and canola. A previous study by Cavalaris et al. 2021 reports similar findings where a model involving seasonal maximum Sentinel-2 EVI had the best performance at prediction of wheat yields. Additionally, our findings on the sensor-specific optimal time period to source predictors, that is at the height of July with Landsat indices and during June-to-August for crop-specific Sentinel-2 indices, could be useful toward developing within-season yield forecasting models. We further note that overall performance of Landsat-based indices were found to be better than MODIS-EVI2 during the 2009–19 period. Such an insight could be valuable for agricultural risk planning and management. Agricultural stakeholders such as crop insurers currently using MODIS-based data for applications such as designing index-based insurance could consider Landsat-based indices to identify any difference in value proposition for farmers when using different satellite data inputs for insurance design. In the near-recent modeling period (2016–19), MODIS-EVI2 had a better prediction performance for canola over all other indices. Such a finding indicates the importance of both better spatial resolution of Landsat and better temporal resolution of MODIS in tracking crop yields at large scales of aggregation. A previous study reported similar findings when MODIS- and Landsat-based indices were tested for monitoring crop yields in contrasting agricultural regions of the U.S. and India (Azzari et al. 2017).

Although crop yield models using Sentinel-2 indices were unable to demonstrate marked improvements over MODIS- and Landsat-based models at the coarser municipality level, many studies have highlighted the potential of Sentinel-2 data for yield prediction at finer spatial levels (e.g., field, farm, landscape) (Hunt et al. 2019; Segarra et al. 2022; Zhao et al. 2020; Cavalaris et al. 2021; Kayad et al. 2019). Our future work will look into evaluation of the best performing Sentinel-2 indices identified in the current study using high resolution field level crop yield data in the Canadian Prairies. For example, we found that green index (GI) outperformed NDWI in tracking plant water signals for wheat and canola in each multi-index combination tested using Sentinel-2 data. GI has been noted to be sensitive to irrigation status in a previous study (Deines et al. 2017) and further testing at the field level could help guide management decisions for yield optimization.

The spatial pattern of yields estimated using different models varied slightly for each crop across some agricultural regions of the Canadian Prairies. These results indicate that chosen satellite index/indices drive the location-specific performance of models requiring further testing and identification of location-specific best vegetation index for monitoring crop health at the sub-regional or field level in the Prairies. For example, Bolton and Friedl 2013 recommends a mixed model, NDWI in more arid regions (for its ability to detect irrigation) and EVI2 in other regions, for predicting maize yield in U.S. counties.

The results of this study showed that statistical yield models based solely on vegetation indices were unable to capture well the extreme training data points. Additionally, model performance was weak in some years shown by the upward and downward spikes in model residuals (Figures A1 and A2). Models highly underestimated crop yields in years 2009 and 2013 and overestimated yields in 2012. A study investigating climate variability impacts on crop production in the Canadian Prairies over the 2000–13 period found that drought conditions in 2009 and 2012 resulted in below average crop productivity, mostly in Alberta and Saskatchewan (Dong et al. 2016). Incorporating additional predictors related to weather and environmental data could help improve performance of models driven solely by vegetation indices in such scenarios. As pure vegetation index-based methods are unable to track yield outcomes in atypical growth conditions (e.g., droughts, high temperatures), a study tested several input datasets (vegetation indices, weather data, environmental variables) and models (linear regression, polynomial regression and machine learning) for estimating yields (Sakamoto 2020). They found that a machine learning model, Random Forest, incorporating additional environmental variables provided accurate yield estimates during a drought year. Statistical models based on machine learning algorithms has advantages over a traditional regression model especially when some predictors are highly correlated (Jeong et al. 2016). Efforts are currently underway to integrate weather and environmental data in satellite-based models of crop yield. Introducing additional information representing variability of environmental conditions may improve the current performance of models in atypical production years. We plan to further evaluate the potential of state-of-the-art machine and deep learning algorithms in detecting spatiotemporal variation in yields over traditional regression models using high spatial resolution yield datasets. Such efforts will provide further exploration of predictors that improve representation of atypical impacts on final crop yield which could be useful for stakeholders such as crop insurers, commodity traders and policy analysts.

Conclusion

Timely and consistent monitoring of crop production using a satellite-based approach is one way to enhance food production and overall food security. A satellite-based approach is advantageous as it can be applied in data sparse regions to obtain crop yield estimates. This study assessed the performance of several vegetation indices derived from different satellite datasets within two modeling approaches for crop yield estimation in the Canadian Prairies. Both single-index and multi-index models demonstrated similar performance at yield estimation over a 11-year period (2009–2019) for each crop (wheat, barley, canola). Single-index, NDWI, outperformed other Landsat-based indices for each crop during this period. Our analysis further showed that best performing crop-specific model using Landsat data produced yield estimates that were marginally better than models using MODIS-EVI2. Spatial examination of estimated yields revealed slightly varying location-specific performance of models trained using different index predictors. Such findings indicate the need to identify the best vegetation index/indices within different agricultural subregions toward improving statistical-based yield estimates. Furthermore, introduction of additional weather and environmental variables within pure vegetation index-based yield models may help improve prediction of interannual yield variability within a given administrative unit or in atypical production years (e.g., droughts, high temperatures). Overall, this study highlights the importance of exploring and testing different vegetation indices capturing different aspects of plant growth toward improving yield estimates and for further use in applications such as yield forecasting in agricultural regions under changing growing conditions.

Acknowledgements

We thank Agriculture Financial Services Corporation (AFSC) for providing crop yield data for Alberta (2009–2019).

Disclosure statement

The authors report that there are no competing interests to declare.

Additional information

Funding

This work was funded by Alberta Innovates Project #2019F163R, entitled “Climate-smart agriculture: Big-data analysis to investigate what management practices work”. NKN was supported by funding from the Canadian Agricultural Partnership (CAP).

Notes

1 AVHRR: Advanced Very High Resolution Radiometer

2 MODIS: Moderate Resolution Imaging Spectroradiometer

3 Saskatchewan crop yields: https://applications.saskatchewan.ca/agrrmyields

4 Manitoba crop yields: https://www.masc.mb.ca/masc.nsf/mmpp_index.html

5 As including multi-step information from the same growing season introduces collinearity between predictors and influences sensitivity of model coefficients, we decided to eliminate multi-stage models from the linear regression modelling approach.

6 Regression models with municipality as a fixed effect were also tested (results not shown). As models with fixed effects did not influence results significantly and due to flexibility of models including a random effect, we present and discuss results for latter models.

Appendix A

 

Figure A1. Accuracy of yearly yield predictions by the best performing multi-index and single-index models using Landsat data against MODIS single-index model for 2009–2019 for wheat, barley and canola. Yield anomaly (%) is the model residual relative to mean observed yield. A simple average of yield anomaly across all municipalities was calculated to derive an annual mean anomaly value.

Figure A2. Accuracy of yearly yield predictions by the best performing multi-index and single-index models using Landsat and Sentinel-2 data against MODIS single-index model for 2016–2019 for wheat, barley and canola. Yield anomaly (%) is the model residual relative to mean observed yield. A simple average of yield anomaly across all municipalities was calculated to derive an annual mean anomaly value.

Figure A3. Correlations between different Landsat-based vegetation indices for wheat observed in July. Pairwise correlation coefficient and pairwise relationship between variables is shown in the upper triangle portion and lower triangle portion of the plot respectively. The curves on the diagonal show the density of respective variables where frequency count is reported along the y-axis and the variable value range is reported along the x-axis.

Figure A4. Correlations between different Landsat-based vegetation indices for barley observed in July. Pairwise correlation coefficient and pairwise relationship between variables is shown in the upper triangle portion and lower triangle portion of the plot respectively. The curves on the diagonal show the density of respective variables where frequency count is reported along the y-axis and the variable value range is reported along the x-axis.

Figure A5. Correlations between different Landsat-based vegetation indices for canola observed in July. Pairwise correlation coefficient and pairwise relationship between variables is shown in the upper triangle portion and lower triangle portion of the plot respectively. The curves on the diagonal show the density of respective variables where frequency count is reported along the y-axis and the variable value range is reported along the x-axis.

Figure A6. Yield estimation accuracy (RMSE) of single stage and multi-stage models using Sentinel-2 data for 2016-2019 for wheat, barley and canola. In the plots, x-axis show a week or combination of weeks using Sentinel-2 data. Single weeks are denoted by labels starting with ‘1’. The week number shown corresponds to a week in a calendar year.

Figure A7. Spatial distribution of observed yield of wheat and estimation error of best-performing models using different satellite index predictors for 2009 and 2019. Landsat-NDWI is the best performing single-index. Best performing multi-index combination using Landsat data for wheat is EVI + SR + NDWI. Regions in white represent the absence of cropland within a given municipality. Regions in grey represent municipalities for which model-specific yield estimation error is unavailable.

Figure A8. Spatial distribution of observed yield of barley and estimation error of best-performing models using different satellite index predictors for 2009 and 2019. Landsat-NDWI is the best performing single-index. Best performing multi-index combination using Landsat data for wheat is EVI + SR + GI. Regions in white represent the absence of cropland within a given municipality. Regions in grey represent municipalities for which model-specific yield estimation error is unavailable.

Figure A9. Spatial distribution of observed yield of canola and estimation error of best-performing models using different satellite index predictors for 2009 and 2019. Landsat-NDWI is the best performing single-index. Best performing multi-index combination using Landsat data for wheat is EVI + SR + NDWI. Regions in white represent the absence of cropland within a given municipality. Regions in grey represent municipalities for which model-specific yield estimation error is unavailable.

Table A1. Cross validation performance metrics for wheat of best performing multi-index models and single indices during the time period corresponding to Sentinel-2 data availability (2016–2019).

	M	L	L	L	L	S	S	S	S
	EVI2	NDVI	EVI	SI	MI	NDVI	EVI	TCARI	MI
RMSE	0.582	0.440	0.441	0.426	0.434	0.651	0.605	0.650	0.711
RRMSE	17.49	13.21	13.24	12.79	13.05	19.56	18.17	19.51	21.35
MAE	0.450	0.341	0.336	0.334	0.339	0.517	0.482	0.513	0.577
$R^2$	0.66	0.64	0.60	0.59	0.57	0.38	0.39	0.37	0.34
WD	0.180	0.129	0.113	0.116	0.120	0.238	0.228	0.283	0.314

The best performing model (in bold) across model types was chosen using RMSE (RRMSE) and corresponding $R^2$ and WD values are bolded for the same model.

M: MODIS-based model; L: Landsat-based model; S: Sentinel-2-based model; SI: best performing single-index model; MI: best performing multi-index model.

SI (L) model: GI.

SI (S) model: EVI.

MI (L) model: EVI + SR + GI.

MI (S) model: EVI + TCARI + GI.

Root-mean-square-error (RMSE) in t/ha, RMSE relative to mean observed yield (RRMSE (%)), mean absolute error (MAE) in t/ha, coefficient of determination (R²), and Wasserstein Distance (WD).

Table A2. Cross validation performance metrics for barley of best performing multi-index models and single indices during the time period corresponding to Sentinel-2 data availability (2016–2019).

	M	L	L	L	L	S	S	S	S
	EVI2	NDVI	EVI	SI	MI	NDVI	EVI	SR	MI
RMSE	0.743	0.719	0.752	0.677	0.718	0.865	0.902	0.851	0.837
RRMSE	16.47	15.93	16.68	15.02	15.91	19.19	20	18.88	18.57
MAE	0.568	0.563	0.596	0.524	0.551	0.698	0.741	0.688	0.685
$R^2$	0.56	0.47	0.44	0.49	0.42	0.26	0.28	0.28	0.29
WD	0.222	0.245	0.269	0.215	0.258	0.436	0.430	0.439	0.456

The best performing model (in bold) across model types was chosen using RMSE (RRMSE) and corresponding $R^2$ and WD values are bolded for the same model.

M: MODIS-based model; L: Landsat-based model; S: Sentinel-2-based model; SI: best performing single-index model; MI: best performing multi-index model.

SI (L) model: NDWI.

SI (S) model: SR.

MI (L) model: EVI + SR + NDWI.

MI (S) model: EVI + SR + NDWI.

Root-mean-square-error (RMSE) in t/ha, RMSE relative to mean observed yield (RRMSE (%)), mean absolute error (MAE) in t/ha, coefficient of determination (R²), and Wasserstein Distance (WD).

Table A3. Cross validation performance metrics for canola of best performing multi-index models and single indices during the time period corresponding to Sentinel-2 data availability (2016–2019).

	M	L	L	L	L	S	S	S	S
	EVI2	NDVI	EVI	NDWI	MI	NDVI	EVI	TCARI	MI
RMSE	0.359	0.435	0.442	0.401	0.416	0.483	0.405	0.487	0402
RRMSE	13.52	16.43	16.67	15.12	15.69	18.23	15.31	18.39	15.19
MAE	0.281	0.343	0.348	0.313	0.326	0.383	0.328	0.381	0.317
$R^2$	0.53	0.35	0.32	0.41	0.37	0.34	0.37	0.28	0.39
WD	0.118	0.168	0.180	0.154	0.162	0.152	0.133	0.100	0.127

The best performing model (in bold) across model types was chosen using RMSE (RRMSE) and corresponding $R^2$ and WD values are bolded for the same model.

M: MODIS-based model; L: Landsat-based model; S: Sentinel-2-based model; SI: best performing single-index model; MI: best performing multi-index model.

SI (L) model: NDWI.

SI (S) model: EVI.

MI (L) model: EVI + SR + NDWI.

MI (S) model: EVI + TCARI + GI.

Root-mean-square-error (RMSE) in t/ha, RMSE relative to mean observed yield (RRMSE (%)), mean absolute error (MAE) in t/ha, coefficient of determination (R²), and Wasserstein Distance (WD).

Table A4. AIC and AICc values of Landsat-based multi-index models for barley, canola, and wheat for 2009–2019 and 2016–2019. The best performing model is presented in bold.

Month	Index combination	Barley 2009–2019	Canola 2009–2019	Wheat 2009–2019	Barley 2016–2019	Canola 2016–2019	Wheat 2016–2019
6	NDVI_GCVI_NDWI	10,211.98	5156.98	6680.69	3554.97	1888.58	2510.60
7	NDVI_GCVI_NDWI	8358.40	3994.57	4707.70	2718.16	1475.02	1580.93
8	NDVI_GCVI_NDWI	9913.34	3898.32	6049.85	3489.41	1340.89	2430.28
6	NDVI_GCVI_GI	10,219.63	5157.93	6778.80	3560.11	1901.23	2549.86
7	NDVI_GCVI_GI	8786.79	4536.27	5098.80	2843.00	1734.61	1612.51
8	NDVI_GCVI_GI	10,308.78	4524.51	6447.14	3556.90	1506.29	2494.40
6	NDVI_GNDVI_NDWI	10,454.46	5151.28	6671.92	3693.07	1908.43	2574.24
7	NDVI_GNDVI_NDWI	9710.15	4170.30	5531.50	3320.82	1521.35	2034.62
8	NDVI_GNDVI_NDWI	10,248.78	4076.91	6126.90	3590.78	1387.44	2422.06
6	NDVI_GNDVI_GI	10,296.12	5088.90	6652.97	3597.38	1896.54	2539.55
7	NDVI_GNDVI_GI	8745.48	4486.66	4970.99	2859.81	1736.55	1614.49
8	NDVI_GNDVI_GI	10,425.78	4396.07	6444.78	3599.46	1506.89	2464.41
6	NDVI_SR_NDWI	10,186.79	5160.45	6466.15	3496.31	1893.36	2401.33
7	NDVI_SR_NDWI	7881.42	3764.34	4475.89	2402.61	1375.17	1535.42
8	NDVI_SR_NDWI	10,019.80	3597.76	6065.65	3518.45	1250.59	2430.35
6	NDVI_SR_GI	10,122.36	5056.05	6456.91	3504.03	1878.09	2434.46
7	NDVI_SR_GI	8033.66	4275.92	4673.52	2478.74	1588.93	1473.35
8	NDVI_SR_GI	10,339.18	4165.83	6403.20	3570.53	1393.33	2464.81
6	EVI_GCVI_NDWI	10,138.14	4642.44	6572.72	3534.26	1733.18	2456.52
7	EVI_GCVI_NDWI	8258.70	3575.84	4532.68	2664.14	1349.13	1522.25
8	EVI_GCVI_NDWI	10,065.31	3686.72	6093.40	3483.45	1220.63	2425.21
6	EVI_GCVI_GI	10,157.00	4651.21	6675.26	3532.58	1732.79	2481.38
7	EVI_GCVI_GI	8433.52	3957.98	4729.63	2795.45	1536.18	1557.50
8	EVI_GCVI_GI	10,246.42	3957.15	6387.11	3537.76	1339.73	2502.78
6	EVI_GNDVI_NDWI	10,357.36	4679.46	6606.14	3631.90	1737.80	2477.31
7	EVI_GNDVI_NDWI	9306.24	3679.71	5161.44	3177.15	1378.59	1877.67
8	EVI_GNDVI_NDWI	10,322.32	3749.96	6165.12	3586.14	1255.72	2424.58
6	EVI_GNDVI_GI	10,260.78	4642.03	6647.56	3566.24	1729.19	2483.04
7	EVI_GNDVI_GI	8398.01	3953.79	4659.00	2806.50	1532.86	1556.30
8	EVI_GNDVI_GI	10,340.46	3871.91	6385.89	3575.55	1340.58	2464.04
6	EVI_SR_NDWI	10,140.60	4658.91	6424.26	3492.09	1717.52	2405.29
7	EVI_SR_NDWI	7889.93	3444.29	4376.58	2376.50	1248.63	1508.21
8	EVI_SR_NDWI	10,143.18	3513.31	6077.95	3515.40	1153.69	2425.71
6	EVI_SR_GI	10,094.52	4578.51	6429.86	3489.31	1712.51	2429.94
7	EVI_SR_GI	7839.94	3687.47	4427.49	2465.04	1404.41	1451.21
8	EVI_SR_GI	10,284.95	3713.22	6352.68	3551.97	1262.39	2475.23

Table A5. AICc values of Sentinel-2-based multi-index models for barley and wheat for 2016–2019. The best performing model is presented in bold.

Month	Week number	Index combination	Barley	Wheat
6	25	NDVI_GCVI_NDWI	1375.61	1156.04
6	26	NDVI_GCVI_NDWI	1666.43	1305.59
7	27	NDVI_GCVI_NDWI	1808.51	1325.89
7	28	NDVI_GCVI_NDWI	1645.04	1266.08
7	29	NDVI_GCVI_NDWI	2098.06	1539.47
7	30	NDVI_GCVI_NDWI	1888.81	1367.50
7	31	NDVI_GCVI_NDWI	2488.36	1611.40
6	25	NDVI_GCVI_GI	1381.51	1158.95
6	26	NDVI_GCVI_GI	1688.40	1376.33
7	27	NDVI_GCVI_GI	1825.33	1376.02
7	28	NDVI_GCVI_GI	1756.86	1489.37
7	29	NDVI_GCVI_GI	2188.19	1789.22
7	30	NDVI_GCVI_GI	2031.73	1577.87
7	31	NDVI_GCVI_GI	2618.30	1866.00
6	25	NDVI_GNDVI_NDWI	1399.77	1163.31
6	26	NDVI_GNDVI_NDWI	1849.23	1368.52
7	27	NDVI_GNDVI_NDWI	1948.30	1443.70
7	28	NDVI_GNDVI_NDWI	1900.08	1383.29
7	29	NDVI_GNDVI_NDWI	2266.32	1658.77
7	30	NDVI_GNDVI_NDWI	2029.89	1488.75
7	31	NDVI_GNDVI_NDWI	2579.37	1684.17
6	25	NDVI_GNDVI_GI	1381.18	1155.77
6	26	NDVI_GNDVI_GI	1649.46	1332.23
7	27	NDVI_GNDVI_GI	1762.15	1321.71
7	28	NDVI_GNDVI_GI	1722.04	1403.90
7	29	NDVI_GNDVI_GI	2177.63	1778.34
7	30	NDVI_GNDVI_GI	2026.84	1563.80
7	31	NDVI_GNDVI_GI	2612.86	1863.44
6	25	NDVI_SR_NDWI	1363.50	1145.64
6	26	NDVI_SR_NDWI	1717.57	1270.52
7	27	NDVI_SR_NDWI	1782.27	1245.66
7	28	NDVI_SR_NDWI	1638.11	1250.45
7	29	NDVI_SR_NDWI	2095.44	1560.15
7	30	NDVI_SR_NDWI	1956.83	1441.94
7	31	NDVI_SR_NDWI	2534.00	1628.66
6	25	NDVI_SR_GI	1372.44	1143.92
6	26	NDVI_SR_GI	1644.97	1300.66
7	27	NDVI_SR_GI	1757.73	1245.77
7	28	NDVI_SR_GI	1682.08	1365.26
7	29	NDVI_SR_GI	2145.92	1759.40
7	30	NDVI_SR_GI	2028.30	1544.21
7	31	NDVI_SR_GI	2614.78	1847.41
6	25	NDVI_NDRE1_NDWI	1378.25	1159.15
6	26	NDVI_NDRE1_NDWI	1844.46	1377.34
7	27	NDVI_NDRE1_NDWI	1837.84	1355.31
7	28	NDVI_NDRE1_NDWI	1945.20	1408.41
7	29	NDVI_NDRE1_NDWI	2313.82	1682.11
7	30	NDVI_NDRE1_NDWI	2174.07	1556.04
7	31	NDVI_NDRE1_NDWI	2638.06	1752.37
6	25	NDVI_NDRE1_GI	1357.85	1146.13
6	26	NDVI_NDRE1_GI	1648.66	1334.24
7	27	NDVI_NDRE1_GI	1721.15	1259.37
7	28	NDVI_NDRE1_GI	1684.09	1359.95
7	29	NDVI_NDRE1_GI	2144.51	1712.61
7	30	NDVI_NDRE1_GI	2017.99	1541.35
7	31	NDVI_NDRE1_GI	2618.22	1840.64
6	25	NDVI_NDRE2_NDWI	1380.79	1165.29
6	26	NDVI_NDRE2_NDWI	1846.19	1381.88
7	27	NDVI_NDRE2_NDWI	1811.26	1327.71
7	28	NDVI_NDRE2_NDWI	1938.04	1401.64
7	29	NDVI_NDRE2_NDWI	2308.75	1679.24
7	30	NDVI_NDRE2_NDWI	2187.11	1548.84
7	31	NDVI_NDRE2_NDWI	2633.82	1753.43
6	25	NDVI_NDRE2_GI	1359.01	1154.54
6	26	NDVI_NDRE2_GI	1658.25	1346.25
7	27	NDVI_NDRE2_GI	1701.59	1256.70
7	28	NDVI_NDRE2_GI	1673.56	1319.67
7	29	NDVI_NDRE2_GI	2135.96	1697.56
7	30	NDVI_NDRE2_GI	2021.75	1520.21
7	31	NDVI_NDRE2_GI	2617.97	1833.59
6	25	NDVI_TCARI_NDWI	1359.47	1096.71
6	26	NDVI_TCARI_NDWI	1882.57	1358.17
7	27	NDVI_TCARI_NDWI	1958.24	1417.69
7	28	NDVI_TCARI_NDWI	2017.61	1432.38
7	29	NDVI_TCARI_NDWI	2252.05	1661.22
7	30	NDVI_TCARI_NDWI	2069.60	1499.94
7	31	NDVI_TCARI_NDWI	2477.19	1675.82
6	25	NDVI_TCARI_GI	1336.17	1089.15
6	26	NDVI_TCARI_GI	1681.42	1338.84
7	27	NDVI_TCARI_GI	1781.00	1336.51
7	28	NDVI_TCARI_GI	1743.75	1478.15
7	29	NDVI_TCARI_GI	2085.15	1750.48
7	30	NDVI_TCARI_GI	1918.34	1539.66
7	31	NDVI_TCARI_GI	2451.71	1797.59
6	25	EVI_GCVI_NDWI	1180.97	955.04
6	26	EVI_GCVI_NDWI	1449.06	1135.16
7	27	EVI_GCVI_NDWI	1544.48	1051.60
7	28	EVI_GCVI_NDWI	1388.59	1003.89
7	29	EVI_GCVI_NDWI	1801.84	1288.58
7	30	EVI_GCVI_NDWI	1796.80	1256.88
7	31	EVI_GCVI_NDWI	2380.38	1458.40
6	25	EVI_GCVI_GI	1192.06	955.11
6	26	EVI_GCVI_GI	1457.00	1162.86
7	27	EVI_GCVI_GI	1544.50	1051.73
7	28	EVI_GCVI_GI	1409.38	1070.91
7	29	EVI_GCVI_GI	1851.17	1327.90
7	30	EVI_GCVI_GI	1873.27	1391.76
7	31	EVI_GCVI_GI	2482.23	1647.39
6	25	EVI_GNDVI_NDWI	1206.21	961.89
6	26	EVI_GNDVI_NDWI	1622.56	1167.85
7	27	EVI_GNDVI_NDWI	1658.55	1156.47
7	28	EVI_GNDVI_NDWI	1576.96	1082.08
7	29	EVI_GNDVI_NDWI	1949.19	1359.72
7	30	EVI_GNDVI_NDWI	1922.13	1379.93
7	31	EVI_GNDVI_NDWI	2454.74	1502.43
6	25	EVI_GNDVI_GI	1190.81	954.59
6	26	EVI_GNDVI_GI	1447.75	1122.21
7	27	EVI_GNDVI_GI	1505.36	1006.47
7	28	EVI_GNDVI_GI	1393.38	1029.98
7	29	EVI_GNDVI_GI	1850.34	1325.16
7	30	EVI_GNDVI_GI	1868.44	1385.77
7	31	EVI_GNDVI_GI	2477.13	1647.39
6	25	EVI_SR_NDWI	1173.68	947.88
6	26	EVI_SR_NDWI	1514.53	1118.79
7	27	EVI_SR_NDWI	1541.51	1023.33
7	28	EVI_SR_NDWI	1385.56	1007.43
7	29	EVI_SR_NDWI	1812.25	1304.35
7	30	EVI_SR_NDWI	1858.61	1332.07
7	31	EVI_SR_NDWI	2409.85	1474.09
6	25	EVI_SR_GI	1189.43	948.14
6	26	EVI_SR_GI	1439.65	1135.04
7	27	EVI_SR_GI	1510.07	1023.34
7	28	EVI_SR_GI	1370.87	1032.11
7	29	EVI_SR_GI	1821.39	1318.75
7	30	EVI_SR_GI	1873.11	1379.96
7	31	EVI_SR_GI	2479.43	1640.31
6	25	EVI_NDRE1_NDWI	1208.74	963.99
6	26	EVI_NDRE1_NDWI	1648.33	1214.23
7	27	EVI_NDRE1_NDWI	1688.80	1177.84
7	28	EVI_NDRE1_NDWI	1682.92	1142.63
7	29	EVI_NDRE1_NDWI	2053.49	1409.57
7	30	EVI_NDRE1_NDWI	2168.17	1482.01
7	31	EVI_NDRE1_NDWI	2539.67	1603.00
6	25	EVI_NDRE1_GI	1186.69	954.58
6	26	EVI_NDRE1_GI	1451.25	1149.77
7	27	EVI_NDRE1_GI	1518.12	1006.78
7	28	EVI_NDRE1_GI	1408.45	1068.99
7	29	EVI_NDRE1_GI	1832.14	1314.02
7	30	EVI_NDRE1_GI	1877.67	1385.49
7	31	EVI_NDRE1_GI	2484.33	1640.87
6	25	EVI_NDRE2_NDWI	1223.72	969.76
6	26	EVI_NDRE2_NDWI	1669.95	1227.57
7	27	EVI_NDRE2_NDWI	1702.44	1198.80
7	28	EVI_NDRE2_NDWI	1740.01	1189.64
7	29	EVI_NDRE2_NDWI	2073.82	1434.75
7	30	EVI_NDRE2_NDWI	2218.09	1509.75
7	31	EVI_NDRE2_NDWI	2544.76	1617.34
6	25	EVI_NDRE2_GI	1185.07	955.11
6	26	EVI_NDRE2_GI	1451.81	1145.77
7	27	EVI_NDRE2_GI	1491.39	993.70
7	28	EVI_NDRE2_GI	1399.28	1067.68
7	29	EVI_NDRE2_GI	1821.11	1310.44
7	30	EVI_NDRE2_GI	1875.96	1378.86
7	31	EVI_NDRE2_GI	2484.62	1638.28
6	25	EVI_TCARI_NDWI	1218.31	945.65
6	26	EVI_TCARI_NDWI	1725.83	1226.91
7	27	EVI_TCARI_NDWI	1774.14	1182.00
7	28	EVI_TCARI_NDWI	1934.52	1286.14
7	29	EVI_TCARI_NDWI	2187.28	1503.04
7	30	EVI_TCARI_NDWI	2149.50	1530.45
7	31	EVI_TCARI_NDWI	2398.16	1588.69
6	25	EVI_TCARI_GI	1180.64	932.35
6	26	EVI_TCARI_GI	1456.68	1144.02
7	27	EVI_TCARI_GI	1523.32	1034.01
7	28	EVI_TCARI_GI	1402.06	1055.24
7	29	EVI_TCARI_GI	1791.59	1306.62
7	30	EVI_TCARI_GI	1780.31	1370.85
7	31	EVI_TCARI_GI	2341.02	1609.56

Table A6. AICc values of Sentinel-2-based multi-index models for canola for 2016–2019. The best performing model is presented in bold.

Month	Week number	Index combination	Canola
6	26	NDVI_GCVI_NDWI	1062.97
7	27	NDVI_GCVI_NDWI	1052.23
7	28	NDVI_GCVI_NDWI	1134.85
7	29	NDVI_GCVI_NDWI	1228.58
7	30	NDVI_GCVI_NDWI	889.56
7	31	NDVI_GCVI_NDWI	1036.94
8	32	NDVI_GCVI_NDWI	834.21
8	33	NDVI_GCVI_NDWI	906.71
8	34	NDVI_GCVI_NDWI	1037.25
6	26	NDVI_GCVI_GI	1059.71
7	27	NDVI_GCVI_GI	1060.42
7	28	NDVI_GCVI_GI	1177.18
7	29	NDVI_GCVI_GI	1426.40
7	30	NDVI_GCVI_GI	1161.71
7	31	NDVI_GCVI_GI	1307.86
8	32	NDVI_GCVI_GI	1009.59
8	33	NDVI_GCVI_GI	1003.53
8	34	NDVI_GCVI_GI	1109.44
6	26	NDVI_GNDVI_NDWI	1111.92
7	27	NDVI_GNDVI_NDWI	1131.05
7	28	NDVI_GNDVI_NDWI	1315.73
7	29	NDVI_GNDVI_NDWI	1349.50
7	30	NDVI_GNDVI_NDWI	1117.07
7	31	NDVI_GNDVI_NDWI	1263.26
8	32	NDVI_GNDVI_NDWI	953.40
8	33	NDVI_GNDVI_NDWI	986.44
8	34	NDVI_GNDVI_NDWI	1106.89
6	26	NDVI_GNDVI_GI	1034.87
7	27	NDVI_GNDVI_GI	1002.87
7	28	NDVI_GNDVI_GI	1157.87
7	29	NDVI_GNDVI_GI	1387.53
7	30	NDVI_GNDVI_GI	1037.78
7	31	NDVI_GNDVI_GI	1186.85
8	32	NDVI_GNDVI_GI	993.96
8	33	NDVI_GNDVI_GI	950.83
8	34	NDVI_GNDVI_GI	1072.32
6	26	NDVI_SR_NDWI	1062.38
7	27	NDVI_SR_NDWI	1039.71
7	28	NDVI_SR_NDWI	1182.60
7	29	NDVI_SR_NDWI	1252.49
7	30	NDVI_SR_NDWI	852.80
7	31	NDVI_SR_NDWI	938.16
8	32	NDVI_SR_NDWI	824.61
8	33	NDVI_SR_NDWI	895.69
8	34	NDVI_SR_NDWI	1041.87
6	26	NDVI_SR_GI	1033.31
7	27	NDVI_SR_GI	1033.52
7	28	NDVI_SR_GI	1156.63
7	29	NDVI_SR_GI	1403.45
7	30	NDVI_SR_GI	998.46
7	31	NDVI_SR_GI	1093.45
8	32	NDVI_SR_GI	957.09
8	33	NDVI_SR_GI	967.90
8	34	NDVI_SR_GI	1096.34
6	26	NDVI_NDRE1_NDWI	1105.10
7	27	NDVI_NDRE1_NDWI	1084.19
7	28	NDVI_NDRE1_NDWI	1339.64
7	29	NDVI_NDRE1_NDWI	1356.02
7	30	NDVI_NDRE1_NDWI	1099.59
7	31	NDVI_NDRE1_NDWI	1287.59
8	32	NDVI_NDRE1_NDWI	947.53
8	33	NDVI_NDRE1_NDWI	1017.87
8	34	NDVI_NDRE1_NDWI	1131.33
6	26	NDVI_NDRE1_GI	1056.72
7	27	NDVI_NDRE1_GI	1055.86
7	28	NDVI_NDRE1_GI	1153.37
7	29	NDVI_NDRE1_GI	1358.93
7	30	NDVI_NDRE1_GI	993.65
7	31	NDVI_NDRE1_GI	1267.10
8	32	NDVI_NDRE1_GI	919.13
8	33	NDVI_NDRE1_GI	975.24
8	34	NDVI_NDRE1_GI	1110.48
6	26	NDVI_NDRE2_NDWI	1093.97
7	27	NDVI_NDRE2_NDWI	1070.02
7	28	NDVI_NDRE2_NDWI	1328.79
7	29	NDVI_NDRE2_NDWI	1358.29
7	30	NDVI_NDRE2_NDWI	1114.50
7	31	NDVI_NDRE2_NDWI	1292.62
8	32	NDVI_NDRE2_NDWI	938.83
8	33	NDVI_NDRE2_NDWI	1018.92
8	34	NDVI_NDRE2_NDWI	1131.14
6	26	NDVI_NDRE2_GI	1054.25
7	27	NDVI_NDRE2_GI	1051.53
7	28	NDVI_NDRE2_GI	1164.03
7	29	NDVI_NDRE2_GI	1352.62
7	30	NDVI_NDRE2_GI	988.90
7	31	NDVI_NDRE2_GI	1268.59
8	32	NDVI_NDRE2_GI	895.95
8	33	NDVI_NDRE2_GI	978.58
8	34	NDVI_NDRE2_GI	1111.04
6	26	NDVI_TCARI_NDWI	1092.53
7	27	NDVI_TCARI_NDWI	1084.84
7	28	NDVI_TCARI_NDWI	1307.37
7	29	NDVI_TCARI_NDWI	1352.60
7	30	NDVI_TCARI_NDWI	1130.49
7	31	NDVI_TCARI_NDWI	1223.94
8	32	NDVI_TCARI_NDWI	916.54
8	33	NDVI_TCARI_NDWI	1007.58
8	34	NDVI_TCARI_NDWI	1123.76
6	26	NDVI_TCARI_GI	1025.45
7	27	NDVI_TCARI_GI	998.37
7	28	NDVI_TCARI_GI	1161.56
7	29	NDVI_TCARI_GI	1421.10
7	30	NDVI_TCARI_GI	1152.31
7	31	NDVI_TCARI_GI	1247.26
8	32	NDVI_TCARI_GI	963.43
8	33	NDVI_TCARI_GI	987.04
8	34	NDVI_TCARI_GI	1079.67
6	26	EVI_GCVI_NDWI	956.21
7	27	EVI_GCVI_NDWI	993.61
7	28	EVI_GCVI_NDWI	1078.67
7	29	EVI_GCVI_NDWI	1164.04
7	30	EVI_GCVI_NDWI	820.42
7	31	EVI_GCVI_NDWI	1033.07
8	32	EVI_GCVI_NDWI	791.33
8	33	EVI_GCVI_NDWI	841.36
8	34	EVI_GCVI_NDWI	998.57
6	26	EVI_GCVI_GI	951.85
7	27	EVI_GCVI_GI	993.67
7	28	EVI_GCVI_GI	1078.43
7	29	EVI_GCVI_GI	1213.09
7	30	EVI_GCVI_GI	946.66
7	31	EVI_GCVI_GI	1094.96
8	32	EVI_GCVI_GI	800.76
8	33	EVI_GCVI_GI	869.24
8	34	EVI_GCVI_GI	1034.38
6	26	EVI_GNDVI_NDWI	1006.27
7	27	EVI_GNDVI_NDWI	1078.99
7	28	EVI_GNDVI_NDWI	1199.60
7	29	EVI_GNDVI_NDWI	1226.60
7	30	EVI_GNDVI_NDWI	1014.58
7	31	EVI_GNDVI_NDWI	1154.38
8	32	EVI_GNDVI_NDWI	893.56
8	33	EVI_GNDVI_NDWI	902.92
8	34	EVI_GNDVI_NDWI	1059.61
6	26	EVI_GNDVI_GI	927.54
7	27	EVI_GNDVI_GI	967.99
7	28	EVI_GNDVI_GI	1017.20
7	29	EVI_GNDVI_GI	1146.31
7	30	EVI_GNDVI_GI	886.05
7	31	EVI_GNDVI_GI	942.24
8	32	EVI_GNDVI_GI	785.13
8	33	EVI_GNDVI_GI	814.15
8	34	EVI_GNDVI_GI	1021.72
6	26	EVI_SR_NDWI	964.87
7	27	EVI_SR_NDWI	985.34
7	28	EVI_SR_NDWI	1107.77
7	29	EVI_SR_NDWI	1162.30
7	30	EVI_SR_NDWI	828.52
7	31	EVI_SR_NDWI	1012.73
8	32	EVI_SR_NDWI	781.13
8	33	EVI_SR_NDWI	828.01
8	34	EVI_SR_NDWI	1009.32
6	26	EVI_SR_GI	929.66
7	27	EVI_SR_GI	977.49
7	28	EVI_SR_GI	1071.90
7	29	EVI_SR_GI	1199.11
7	30	EVI_SR_GI	852.00
7	31	EVI_SR_GI	942.46
8	32	EVI_SR_GI	777.77
8	33	EVI_SR_GI	844.46
8	34	EVI_SR_GI	1033.28
6	26	EVI_NDRE1_NDWI	1018.47
7	27	EVI_NDRE1_NDWI	1097.29
7	28	EVI_NDRE1_NDWI	1236.04
7	29	EVI_NDRE1_NDWI	1231.45
7	30	EVI_NDRE1_NDWI	1054.60
7	31	EVI_NDRE1_NDWI	1185.55
8	32	EVI_NDRE1_NDWI	902.36
8	33	EVI_NDRE1_NDWI	937.04
8	34	EVI_NDRE1_NDWI	1082.78
6	26	EVI_NDRE1_GI	940.95
7	27	EVI_NDRE1_GI	979.13
7	28	EVI_NDRE1_GI	1063.75
7	29	EVI_NDRE1_GI	1215.85
7	30	EVI_NDRE1_GI	906.36
7	31	EVI_NDRE1_GI	1092.23
8	32	EVI_NDRE1_GI	802.94
8	33	EVI_NDRE1_GI	866.28
8	34	EVI_NDRE1_GI	1030.25
6	26	EVI_NDRE2_NDWI	1021.17
7	27	EVI_NDRE2_NDWI	1101.20
7	28	EVI_NDRE2_NDWI	1246.84
7	29	EVI_NDRE2_NDWI	1235.34
7	30	EVI_NDRE2_NDWI	1051.73
7	31	EVI_NDRE2_NDWI	1178.93
8	32	EVI_NDRE2_NDWI	901.19
8	33	EVI_NDRE2_NDWI	935.78
8	34	EVI_NDRE2_NDWI	1083.94
6	26	EVI_NDRE2_GI	937.36
7	27	EVI_NDRE2_GI	974.01
7	28	EVI_NDRE2_GI	1057.07
7	29	EVI_NDRE2_GI	1215.76
7	30	EVI_NDRE2_GI	904.34
7	31	EVI_NDRE2_GI	1088.92
8	32	EVI_NDRE2_GI	804.01
8	33	EVI_NDRE2_GI	862.28
8	34	EVI_NDRE2_GI	1030.23
6	26	EVI_TCARI_NDWI	982.62
7	27	EVI_TCARI_NDWI	1031.69
7	28	EVI_TCARI_NDWI	1201.89
7	29	EVI_TCARI_NDWI	1223.71
7	30	EVI_TCARI_NDWI	1023.45
7	31	EVI_TCARI_NDWI	1125.87
8	32	EVI_TCARI_NDWI	869.07
8	33	EVI_TCARI_NDWI	920.42
8	34	EVI_TCARI_NDWI	1074.16
6	26	EVI_TCARI_GI	933.38
7	27	EVI_TCARI_GI	954.18
7	28	EVI_TCARI_GI	1066.26
7	29	EVI_TCARI_GI	1194.76
7	30	EVI_TCARI_GI	927.44
7	31	EVI_TCARI_GI	1025.60
8	32	EVI_TCARI_GI	759.59
8	33	EVI_TCARI_GI	849.38
8	34	EVI_TCARI_GI	1026.59

Table A7. Sensor-specific multi-index combinations tested for modeling.

Sensor	Index combination
Landsat (7 & 8), Sentinel-2	1. NDVI + GCVI + NDWI
	2. NDVI + GCVI + GI
	3. NDVI + GNDVI + NDWI
	4. NDVI + GNDVI + GI
	5. NDVI + SR + NDWI
	6. NDVI + SR + GI
	7. EVI + GCVI + NDWI
	8. EVI + GCVI + GI
	9. EVI + GNDVI + NDWI
	10. EVI + GNDVI + GI
	11. EVI + SR + NDWI
	12. EVI + SR + GI
Sentinel-2	13. NDVI + NDRE1 + NDWI
	14. NDVI + NDRE1 + GI
	15. NDVI + NDRE2 + NDWI
	16. NDVI + NDRE2 + GI
	17. NDVI + TCARI + NDWI
	18. NDVI + TCARI + GI
	19. EVI + NDRE1 + NDWI
	20. EVI + NDRE1 + GI
	21. EVI + NDRE2 + NDWI
	22. EVI + NDRE2 + GI
	23. EVI + TCARI + NDWI
	24. EVI + TCARI + GI

Index combination specification follows, biomass index + chlorophyll index + water index.