Detection and attribution of cereal yield losses using Sentinel-2 and weather data: A case study in South Australia

Weather extremes affect crop production. Remote sensing can help to detect crop damage and estimate lost yield due to weather extremes over large spatial extents. We propose a novel scalable method to predict in-season yield losses at the sub-field level and attribute these to weather extremes. To assess our method ’ s potential, we conducted a proof-of-concept case study on winter cereal paddocks in South Australia using data from 2017 to 2022. To detect crop growth anomalies throughout the growing season, we aligned a two-band Enhanced Vegetation Index (EVI2) time series from Sentinel-2 with thermal time. The deviation between the expected and observed EVI2 time series was defined as the Crop Damage Index (CDI). We assessed the performance of the CDI within specific phenological windows to predict yield loss. Finally, by comparing instances of substantial increase in CDI with different extreme weather indicators, we explored which (combinations of) extreme weather events were likely responsible for the experienced yield reduction. We found that the use of thermal time diminished the temporal deviation of EVI2 time series between years, resulting in the effective construction of typical stress-free crop growth curves. Thermal-time-based EVI2 time series resulted in better prediction of yield reduction than those based on calendar dates. Yield reduction could be predicted before grain-filling (approximately two months before harvest) with an R 2 of 0.83 for wheat and 0.91 for barley. Finally, the combined analysis of CDI curves and extreme weather indices allowed for timely detection of weather-related causes of crop damage, which also captured the spatial variations of crop damage attribution at sub-field level. The proposed framework provides a basis for early warning of crop damage and attributing the damage to weather extremes in near real-time, which should help to adopt appropriate crop protection strategies.


Introduction
Weather extremes, such as drought, frost, flooding, and heat waves, can severely impact agricultural production and threaten food security worldwide (Schmidhuber and Tubiello, 2007;Vogel et al., 2019).Extreme events have been projected to intensify and become more frequent under future climate scenarios (IPCC, 2022;Kharin et al., 2013;Milly et al., 2002).Global warming, shifts in atmospheric circulation, and altered precipitation patterns, contribute to exacerbating the occurrence and severity of these events (Pachauri et al., 2014;Sun et al., 2015).To develop effective strategies for adaptation and mitigation and reduce the adverse effects of weather calamities on agricultural production, it is essential to comprehend the specific responses of crops to extreme weather conditions (Chavez et al., 2015).
Crop production is vulnerable to the capricious impact of weather extremes, with the potential for multiple stressors to occur within a single growing season.The adverse effects of weather extremes can be partially compensated through an adaptation of farmers' practices, such as improved water use efficiency, modifying sowing and harvest schedules, and selecting cultivars that are better suited to the changing conditions (Lobell et al., 2008;Tao et al., 2014).However, optimum adaptation choices hinge upon a comprehensive understanding of the extent and cause of crop losses incurred.Not all stress events that occur within the growing season lead to crop yield loss (Leonard et al., 2014;Zscheischler and Seneviratne, 2017).Crop damage severity can be highly variable within a field, between fields, and across years, emphasizing the necessity for detailed, accurate, and timely insights into crop responses to stressors.Understanding the cause of crop losses is valuable for early warning systems, as well as improving the definition and implementation of resilient agricultural management options to achieve higher or stabler crop yields.
To implement effective early warning systems, understanding when stresses occur is critical, because crop susceptibility to weather extremes changes depending on phenological period (Sánchez et al., 2014;Satorre and Slafer, 1999).For example, if winter wheat flowers too late, it becomes susceptible to drought and heat stress, whereas early flowering increases its vulnerability to frost-induced damage (Zheng et al., 2012).Frost can strongly affect cereal production around anthesis, engendering sterility and the abortion of produced grains (dos Santos et al., 2022).However, frost during the dormant phase can confer benefits to winter cereals, inducing dormancy and facilitating the accumulation of sugars that function as antifreeze agents (Fuller et al., 2007).In parallel, heat stress during the vegetative stage precipitates leaf wilting and diminished photosynthesis, whereas during the reproductive stage, it can reduce pollen viability and disrupt pollination, leading to reduced yields (Stone and Nicolas, 1995;Tashiro and Wardlaw, 1990).Consequently, to unravel the intricate interactions between weather extremes and crop yield, monitoring crop growth anomalies and tracking dynamic changes in crop phenology are important.
Existing research on crop damage assessment relies mostly on crop model simulations with emphasis on the impact of single weather extremes (Barlow et al., 2015;Eitzinger et al., 2013;X. Li et al., 2015;Sutka, 1994).Crop models simulate crop growth under various environmental conditions based on the understanding of genotype, environment, and management interactions (Barnabás et al., 2008).Despite improvements in crop modelling approaches, existing models still fail to adequately capture the impacts of weather extremes on crop growth (Barlow et al., 2015;Rötter et al., 2018).One limitation lies in the simplified response processes in different phenological stages to extreme weather events (Eitzinger et al., 2013;Mittler, 2006).For example, most crop models integrate frost damage via winterkill functions, seedling death, or advanced senescence, while ignoring the disproportionate reduction in the number of grains due to sterility and abortion (Brisson et al., 2003;Jones et al., 2003;Zheng et al., 2014).Moreover, existing research efforts predominantly assess crop damage for single extreme weather influences, with few focusing on the concurrent or consecutive impacts of multiple weather events (Asseng et al., 2011;Barnabás et al., 2008;Challinor et al., 2005).Previous studies have shown that the molecular and metabolic responses of plants to compound stress events are unique from individual events, posing challenges for models to quantify these synergistic effects (Mittler, 2006;Rizhsky et al., 2004).Additionally, uncertainties in model parameterization and the large data demand for calibration purposes contribute to potential inaccuracies in simulated results, especially at large scales.
Remote sensing (RS) can help in understanding yield losses, given its ability to monitor crop growth status for large areas with both high spatial and temporal resolution (Berger et al., 2022;Karthikeyan et al., 2020;Weiss et al., 2020).Vegetation indices (VIs), derived from spectral reflectance measurements, are frequently used to evaluate changes in vegetation abundance and vigour.Multiple studies have documented the close relationship between VIs and crop yields at different spatial scales (Bolton and Friedl, 2013;Johnson, 2014;Kastens et al., 2005;Liu et al., 2018;Shanahan et al., 2001;Wang et al., 2020b;Xiao et al., 2024).Typically, temporal profiles of VIs are compared between the current year (under stress) against historical years (representing normal conditions) (Funk and Budde, 2009;Unganai and Kogan, 1998;Wang et al., 2020a).The resulting anomaly represents a relative vegetation condition, providing a quantitative measure of the damage severity.Given that weather patterns also affect the timing of crop phenological stages, substantial spatial and interannual variability occurs in crop-specific VI time series, irrespective of the occurrence of extreme weather events.This natural variability makes it challenging to define a reference temporal VI profile for a crop that does not experience stress conditions throughout the growing season.Consequently, methods are needed that account for the natural weather variability to disentangle which crop damage-related VI reductions are caused by extreme events.
Despite existing efforts on crop damage assessment, the potential of RS to assist in timely attribution of crop damage to weather extremes remains under-explored.Previous studies on attribution analysis primarily utilize data-driven approaches employing phenologically and spatially dynamic extreme weather indicators (Lüttger and Feike, 2018;Riedesel et al., 2023;Schmitt et al., 2022).These approaches explain yield losses by establishing statistical relationships between extreme weather indicators and end-of-season yield losses.However, due to limitations such as the availability of high-resolution weather data and extensive phenology and yield information required, this method is mainly restricted to post-season assessment at the regional level (Nóia, 2023;Vogel et al., 2019;Zhao et al., 2022).In contrast, high resolution (≤ 10 m and 5 days interval) RS data such as Sentinel-2 time series presents an opportunity for timely detection of crop damage at sub-field level.High temporal and spatial resolution RS imagery, weather data and observed yields, allow us to attribute crop damage to specific extreme weather events by analysing the correlation between the occurrence of observed crop anomalies and stressors.Overall, we still lack a precise understanding of how weather extremes affect crop production, and to what extent the timing of such extremes influences yield reduction.
In this study, we propose a novel and scalable method incorporating both Sentinel-2 time series and meteorological data to predict in-season yield reduction at the sub-field level and attribute the crop damage to different weather extremes.We aim to address the following three research questions: • How can a reference VI profile be constructed that represents a stress-free crop growth curve considering the spatial and temporal variability of the environmental conditions?• How accurately and how early in the season can the RS time series predict within-field reduction of crop yield by capturing stressinduced changes in crop development?• Can the weather-related causes of crop damage be inferred at subfield level within a useful time-frame for growers by combining RS time series with meteorological data?

Study area
Winter cereals, including winter wheat and barley, account for 58% and 22% of the total national grain production in Australia, respectively.We focused on crop paddocks (a uniform description of fenced crop fields in Australia) cultivated with wheat and barley in an area surrounding Jamestown in the Mid North region of South Australia (Fig. 1), which is a highly-productive agricultural region and part of the Australian wheatbelt.It has a Mediterranean climate with cold and wet winter (June to August) during which winter cereals are grown, and hot and dry summer (December to February) during which land is left fallow.Wheat and barley are usually sown in April-May (austral autumn) and mature at the end of spring (November).Fig. 2 shows the different phenological stages of winter cereals in Australia based on the Zadoks decimal growth scale (Stapper, 2007).The Mid North region experiences significant fluctuation in temperature and rainfall within and across growing seasons (as shown in Fig. 3), resulting in substantial variation in crop yields (Asseng et al., 2011).Extreme weather (drought, frost, and heat) is known to regularly cause crop stress in the region, but in recent years stresses have become more frequent and severe (Borchers Arriagada et al., 2020;Collins et al., 2000).

Sentinel-2 data and preprocessing
Temporal profiles of spectral reflectance were extracted from the Sentinel-2 images for six years (2017)(2018)(2019)(2020)(2021)(2022) for the March to December growing season.ESA's Level-2A Bottom-Of-Atmosphere (BOA) products were only available in Google Earth Engine (GEE) Catalogue from December 2018 onwards and consequently did not match the full timeframe.175 Level-2A (L2A) images for tile 54HTJ were obtained from the GEE collection "COPERNICUS/S2_SR" by filtering out images with cloud coverage above 75%.For dates prior to December 2018, we complemented this with 75 Level-1C (L1C) Top-Of-Atmosphere (TOA) images that were downloaded from the Copernicus Open Access Hub.The Sen2Cor processor (version 2.10) was applied for atmospheric correction of each scene to produce equivalent L2A products.For all 2017-2022 images, the Scene Classification Layer (SCL) was applied to remove pixels classified as: no data (0), saturated or defective (1), dark area (2), cloud shadow (3), cloud medium (8) and high ( 9) probability with a 20 m buffer to further reduce the cloud and snow effects.We masked out bright anomalies (e.g., undetected clouds and smoke) following a temporal filter approach adapted from the MAJA cloud detection algorithm (Bolton et al., 2020;Hagolle et al., 2017).This approach captures steep increases in the blue band caused by bright surfaces which could mean missed clouds and smoke.These masked images were then clipped to the study area extent.

Cropland data
Local growers provided annual cropland information from 2017 to  2022, including paddock boundaries (depicted as black polygons in Fig. 1) and crop metrics such as type, rotation, and yield data. i) Crop paddocks: We focused on the paddocks where wheat and barley were cultivated.In these winter cereal paddocks, growers generally use crop rotations such as oil seeds, pulse crops, improved pasture and fallow, or livestock grazing to reduce diseases and weeds and maintain soil health.Our dataset contained 67 wheat paddock years and 68 barley paddock years with sizes ranging from 2 to 80 ha for 2017-2022, where a paddock year represents data obtained for a single paddock for a single year and crop type.
ii) Crop yield data: we obtained 24 yield maps at 10 m resolution for wheat (14 paddock years) and barley (10 paddock years).The yield maps were created from on-the-go data collected from combine harvesters equipped with yield-monitoring sensors.Subsequently, on-thego yield data was calibrated and mapped following the wellestablished yield mapping protocol (Bramley and Jensen, 2014;Taylor et al., 2007) as implemented in the Precision Agriculture Tools (PAT) plugin (Ratcliff et al., 2020) to produce 5 m resolution yield and standard error maps.The yield maps were then reprojected and resampled to 10 m resolution and aligned with the Sentinel-2 image pixels.
iii) Crop calendar: Crop calendar information was collected from grower interviews.From this information, we used sowing dates which were available for five wheat and six barley paddock years.

Meteorological data
Meteorological data were acquired from the SILO (Scientific Information for Land Owners) climate database (https://www.longpaddock.qld.gov.au/silo), which provides gridded weather maps with daily temporal resolution and 0.05 degrees (~5 km) spatial resolution.We used the daily precipitation, minimum, and maximum temperature from 2017 to 2022 from the 12 grid cells that covered the study area.The gridded data were created by spatially interpolating all available weather stations from the Bureau of Meteorology and other providers using ordinary kriging for rainfall and a thin plate smoothing spline for temperature data (Jeffrey et al., 2001).In the study area, despite a modest spatial variation in daily maximum and minimum temperature and precipitation with an average Standard Deviation (SD) of 0.29 • C, 0.18 • C, and 0.63 mm, a substantial temporal disparity was observed (Fig. S2), revealing a notable average SD of 4.14 • C, 3.21 • C, and 2.39 mm between successive years in 2017-2022.

Auxiliary data
Annual crop and pasture reports (CPR) (https://www.pir.sa.gov.au)published by the Government of South Australia provide estimates of the area and production of each crop and summarize the weather conditions and crop performance at the district level.The CPR reports of Spring Crop Performance in the Upper North district together with expert knowledge acquired from local agronomists were applied to provide evidence information about the actual external stresses that have been reported to cause yield reduction of winter cereals in each growing season.

Methodology
We developed the Crop Damage Index (CDI) based on the two-band Enhanced Vegetation Index (EVI2) time series and weather data to assess crop status throughout the growing season.The CDI is an anomaly measure that assesses the difference between the EVI2 time series derived from Sentinel-2 at any given time during the season and a predefined reference EVI2 time series.Considering that climatic conditions and management practices differed between locations and years, a single reference curve derived from the historical time series is not directly comparable with the target curve.Thus, we applied two methods to account for these various factors: i) we utilized accumulated growing degree days (AGDD) to define a thermal time to replace the calendar dates; ii) we used the shape model fitting (SMF) method to adjust the reference curve and match the target curve.Subsequently, we used regressions between the CDI, for different moments in the season, against measured yield (which we transformed into yield reduction).To find the optimal way of calculating CDI, we compared the performance of date and AGDD-based CDI curve, with or without the SMF method on yield reduction assessment.Extreme weather impacts on crops were modelled as moments during the season when CDI increases; we linked such moments to potential stress events that we derived from weather data.The diagram in Fig. 4 presents the workflow that we followed and Steps 1 to 4 are further described in Section 3.1 to 3.4, respectively.

Generating date-based EVI2 time series
EVI2 is designed to enhance the vegetation signal by de-coupling the vegetation from the canopy background signal and reducing atmospheric influences (Huete et al., 1994;Jiang et al., 2008).It is more responsive than the most commonly used normalized difference vegetation index (NDVI) to canopy structural variations and in reducing saturation effects of densely vegetated surfaces.Furthermore, it improves the consistency between different sensors compared with EVI because the measured blue band reflection may vary greatly across various sensors when atmospheric influences are significant; this enhanced consistency may facilitate applying our methods to alternative (or combined) satellite datasets, especially the sensor without blue band.We calculated EVI2 from the pre-processed Sentinel-2 images following Eq.(1) for each pixel within the paddock polygons (Jiang et al., 2008): where Red and NIR refer to the reflectance in the 10 m Sentinel-2 images for spectral bands 4 and 8, respectively.We applied a 20 m inner buffer to remove possible interference with reflection from other land covers around the paddock edges.For each year (i.e., 1 March to 31 December, covering the growing season of wheat and barley), we produced an image stack in which every layer represents a single date of observation.The EVI2 time series were extracted from image stacks for each pixel and resampled to 5-day intervals with linear interpolation.

Constructing thermal-time-based EVI2 time series
Accumulated heat over the growing season is an important factor that affects the crop growth rate and phenological development (Wang, 1960).Heat accumulation is often expressed using the AGDD concept, defined as the accumulation of daily average temperature since the date of sowing or leaf emergence.It has been widely accepted as a key driver of plant phenology and used in crop models to define different crop growth stages from leaf emergence to senescence or harvest (Bonhomme et al., 1994;Lu et al., 2001).It also shows potential in providing better consistency and coherence of temporal VI trajectories in different locations or years when transforming the time scale of RS-derived time series to AGDD (Duveiller et al., 2013).We applied AGDD to replace the time scale in order to construct EVI2 time series that depict crop growth Three sequential steps were carried out to build the time series using thermal time: i. Estimating the onset date of the crop growing season.The onset date can be assessed with VI time series and is a proxy for the date of sowing or leaf emergence.To select the most effective way for retrieving the onset date of the growing season, we tested three commonly used change detection approaches, i.e., the moving average, thresholdbased, and largest derivative method.The analysis and comparison of results are presented in supplementary materials Section S3.Among the three algorithms, the threshold-based method performed best and was used for onset date detection.The threshold is applied following the fitting of a smooth function to the EVI2 time series.For this purpose, the double logistic function was used, as it translates into a monotonic description of green-up and senescence.In this way, we could use a small threshold value to find the moment of initial green-up, which is expected to happen shortly after sowing.
Since the EVI2 has a bias towards lower values due to undetected atmospheric effects, the double logistic function was applied to the upper envelope of the data following an iterative adaption of weights as proposed by Chen et al., 2004.Following the function fit, the onset date is retrieved as the date when the fitted curve value crosses a specific threshold (i.e., a percentage of the amplitude).The threshold was set as 0.5% of the ascending curve, to make it as close as possible to the initial green-up that is close to the sowing date.ii.Calculating AGDD.AGDD is calculated by accumulating the growing degree days (GDD) starting from the onset date detected in the previous step until the end of growing season in December (Eq.2).GDD is estimated as the average of the daily maximum (T max ) and minimum temperatures (T min ) minus a base temperature (T base ) (McMaster & Wilhelm, 1997) (Eq.3).To mitigate the influence of extreme temperature, T base (equal to 0 • C) and an upper temperature limit T cutoff (equal to 25 • C) were applied for both wheat and barley to define the temperature range within which plant growth is considered to progress (Bauer et al., 1990;Cao and Moss, 1989;McMaster and Smika, 1988;McMaster and Wilhelm, 1997): iii.Building thermal-time-based EVI2 time series.The date-based time series were then converted to thermal time by replacing day-ofyear (DOY) with AGDD as the temporal unit of the EVI2 time series.

Establishing a reference EVI2 curve
To quantitatively analyse the EVI2 response to weather extremes, we first constructed two reference curves, one for wheat and one for barley, which can represent the typical growth cycle of these crops under normal, stress-free conditions.As we lacked ground data on whether a specific location experienced no-stress conditions during a particular growing season, we adopted an alternative approach to establish the reference curve.The main assumption is that higher EVI2 values correspond to "no-stress" conditions, when comparing different wheat or barley paddocks in different years.First, EVI2 curves for each pixel within the crop paddocks were extracted and smoothed by applying the Savitzky-Golay (SG) filter with a polynomial order of 2 and a window size of 11 observations iteratively to the upper envelope.To extract the reference curve that represents no-stress conditions, we determined the 80 th to 90 th percentile values of the EVI2 for each calendar date, based on all pixels for which we had observations (i.e., all pixels in the 67 wheat and 68 barley paddock years), and then calculated the average value between the 80 th and 90 th percentile in each date and AGDD unit (Fig. 5a).The interval intends to include most of the pixels from healthy and well-maintained crop areas, mask out low EVI2 values when crops are likely under stress and extremely high EVI2 values that could represent outliers, possibly related to high weed presence.To establish an interval range that can lead to appropriate reference curves, we performed a sensitivity analysis by slightly adapting the percentile values.Besides the average of the 80 th -90 th percentile values, we also tested the average of 85 th -95 th percentile, 75 th -95 th percentile, and 85 th -90 th percentile.As the standard deviations between the produced reference curves were small (range by 0.03-0.08),we limited the analysis to using the average EVI2 values of 80 th -90 th percentile interval.

Crop phenometrics extraction
Using the raw EVI2 time series for all pixels within paddocks (i.e., the observed curves), penalized cubic splines were applied to fit the interpolated date-based and thermal-time-based time series to retrieve metrics that characterize the expression of a phenological period (i.e., phenometric).In contrast to the double logistic function used for globally smoothing curves to detect onset dates, this method captures the inseason anomalies effectively while offering smoother fits than the SG filter by utilizing piecewise polynomial functions, thus enhancing flexibility and adaptability for modelling complex phenological transitions.Because noise and data gaps in time series can negatively influence the correct reconstruction of interpolated daily values (and therefore of the estimated phenometrics), we incorporated data from additional years outside of the target year to reduce the impact of both sources of uncertainty (Bolton et al., 2020).First, a cubic spline was fitted to the EVI2 time series for each pixel.Then we calculated the weights following Eq.4 based on the Euclidean distance between the fitted time series for the target year and each of the other years.
where D y is the Euclidean distance between the fitted EVI2 time series of target year t and each of the other years y, and D max is the maximum Euclidean distance considered.MaxW refers to the maximum weight for the target year and was set as 0.1 following Bolton et al., 2020.The weights were assigned to each of the other years to fit the spline for the target year.
Multiple thresholds were applied to the fitted curves to retrieve phenometrics from EVI2 time series.In total, 11 phenometrics were extracted, named Start of Season (SOS) and End of Season (EOS): SOS10, SOS30, SOS50, SOS70, SOS90, EOS10, EOS30, ESO50, ESO70, EOS90, and PS (Peak Season), which refer to the date when EVI2 first reaches 10%, 30%, 50%, and 90% of the difference between the maximum and minimum EVI2 in green-up (upward) phase and senescence (downward) phase, as well as the date of maximum EVI2 value (Fig. 5b).The retrieved phenometrics were extracted to capture the dynamic changes of phenological periods in relation to crop damage for constructing the CDI curve and better predicting yield reduction.

Shape model fitting
The established reference curve may not perfectly depict the stressfree growth trajectory of the crop in the target paddock or pixel due to variations in climatic conditions and agricultural practices.To address this, we employed the SMF method to adjust the reference curve (as shown in Fig. 5c), aligning it more closely with the EVI2 time series observed in the pixel during the early growth phase, when plants are less susceptible to external pressures − specifically, the early vegetative stage between onset date and SOS30 (Sakamoto, 2018;Sakamoto et al., 2010).The details of SMF method and its implementation can be found in Section S4 of the supplementary materials.The resulting transformed Shape Model (SM) for each wheat and barley pixel was deemed a reasonable approximation of the typical "no-stress" growth curve anticipated at that specific location.

Creating the CDI curve
Fig. 5d displays the SG-filtered EVI2 curve of the target pixel and the transformed reference curve obtained using the SMF method.Extreme weather can affect crop development, which will be manifested as reduced photosynthesis, wilted leaves, sterility, and possibly plant death.If the crop is affected, this will then result in slower growing or declined EVI2 after the stress occurs.The CDI curve is calculated as the difference between the transformed reference curve and the observed curve within the growing season (Fig. 5d).Positive differences were assumed to relate to crop damage.Negative CDI values throughout the season were set to zero.Accordingly, the area under the CDI curve provides a measure of accumulated crop damage throughout the season.

Preparation of input variables
To conduct the sub-field analysis, we divided each paddock-year that had available yield map data (n = 24, 14 for wheat and 10 for barley) into different segments.This division was based on the farming-by-yield approach that is used for delineating homogenous management zones in precision agriculture (Basso et al., 2007;Lark, 1998).This approach uses yield as a proxy for soil variables, terrain differences, and small-scale weather variability to divide the paddock into zones that have different productivity levels (Maestrini and Basso, 2021).The identification of spatially coherent regions was carried out by applying a quantile-based clustering method that classified each paddock year into seven segments based on the yield maps. Accordingly, for each growing season, seven segments were derived per paddock, resulting in a total of 98 segments for wheat and 70 segments for barley for subsequent regression analysis.
With the retrieved phenometrics (Section 3.2.2),we calculated − for each pixel − the CDI value corresponding to each of the 11 phenometrics, as well as the areas under the CDI curve (cumCDI) between each pair of adjacent phenometrics.Calculating cumCDI solely for neighbouring phenometrics, rather than all possible phenometric combinations, aims to diminish variable correlation in regression analysis and intends to capture the impact of stressors on individual short phenological periods.In total, we produced 21 input variables for the regression model, including 11 CDI variables and 10 cumCDI variables.These pixel-level values were then spatially aggregated to the paddock segments by taking the average value of the pixels contained in that segment.

Single variable regression model
To assess the ability of CDI and cumCDI at different phenological stages to predict yield reduction, we first conducted a single variable regression analysis to assess if a single CDI (linked to the timing of one phenometric) or cumCDI can explain yield reduction throughout the study area and analysed years.We applied linear and quadratic regression models to capture both linear and nonlinear relationships.We utilized the 21 phenological metrics derived in Section 3.3.1 as the input variables.Yield reduction was defined as the decrease in crop yield compared to the expected or potential yield under optimal growing conditions.The expected stress-free yield was calculated as the average yield value of pixels which have the 80-90% highest yield in all available yield maps of wheat or barley paddocks (following the same logic as the reference curve construction in Section 3.2.1).To avoid negative yield loss, we set the yield reduction as zero for segments which actual yield above the stress-free yield.For each model, we tested the performance of date-based and AGDD-based EVI2 time series.The coefficient of determination (R 2 ) and normalized root mean squared error (NRMSE) were used as test measures to assess the model performance.

Covariate selection
The previous analysis aimed to find the single most sensitive moment or period for yield reduction estimation, while a combination of CDI for the different phenological stages may potentially predict yield reduction better.Considering the possible multicollinearity among the variables, we used the Pearson correlation coefficient (r) and random forest (RF) feature importance ranking to remove the redundant and irrelevant variables prior to regression analysis.The removal of redundant variables aims to improve model efficiency and reduce the risk of overfitting.
Firstly, we calculated Pearson's r between each pair of the CDI and cumCDI variables to assess possible collinearity.Values of r above 0.8 with statistical significance (p ≤ 0.01) suggest collinearity.The feature importance ranking indicates the relative contribution of each variable to the prediction of yield reduction.We used the RF method considering its ability in managing high-dimensional input variables and its robustness to over-fitting compared with other feature selection algorithms (Pal and Foody, 2010).The final importance of each variable is determined by the permutation importance, which is calculated based on a metric that captures the increase in mean square error (MSE) from out-of-sample predictions after randomly permuting the values of the respective predictor (Goldstein et al., 2011).We utilized an ensemble of 500 decision trees with a maximum depth of 10 levels for the RF regressor to balance model complexity and predictive accuracy.The number of variables considered at each split was set to the square root of the total number of input features as commonly recommended, ensuring diversity among the trees (Belgiu and Drȃgut ¸, 2016).The RF model was repeated 50 times.We divided the dataset (segment-based) into training and testing sets, with a 70:30 ratio.The training data was used to build the RF model, while the test set served for cross-validation.Eventually, the feature importance was calculated as an average of the 50 runs using the Scikit-Learn package in Python 3.7.
According to the ranking order of importance scores, the performance of different numbers of input features for predicting yield reduction was assessed.The retained number of features for further analysis was then selected as the lowest number for which the R 2 curve saturated (i.e., only showed a negligible increase).In this case, the threshold is set to stop the variable number selection at which the increase of R 2 from n to n+1 is less than 5%.Among the selected variables with high importance scores, we identified a subset of non-correlated CDI and cumCDI variables by removing one variable from the pair which shows a high correlation (r ≥ 0.8) with statistical significance (p ≤ 0.01).The removed variables corresponded to the later phenological stage to facilitate early prediction of yield reduction with the retained variables.

Multivariate regression model
Multivariate regression analysis was performed with the RF model using the retained variables from Section 3.3.3as input.RF is an ensemble learning method that combines multiple decision trees to provide robust predictions (Breiman, 2001).It has been used effectively for global and regional scale crop yield predictions (Jeong et al., 2016).
The parameter settings of RF regressor model and the cross-validation assessment are the same as described in Section 3.3.2.To select the optimal CDI construction method, we compared the performance of: i) DOY-EVI2 time series without SMF; ii) AGDD-EVI2 time series without SMF; iii) DOY-EVI2 time series with SMF; and iv) AGDD-EVI2 time series with SMF on yield reduction prediction.
Besides, we also evaluated the robustness of the proposed method to the choice of VI and statistical learning method by comparing 1) the EVI2 time series against three other VIs, i.e., normalized difference vegetation index (NDVI) (Kriegler et al., 1969), normalized difference water index (NDWI) (Gao, 1996), and normalized difference red edge index (NDRE) (Maccioni et al., 2001), each possessing distinct spectral characteristics relevant to crop health status; 2) the RF method against two other widely used regression models, including support vector regression and multiple linear regression, which capture relationships between variables differently than RF (Van Klompenburg et al., 2020).

Extreme weather index calculation
To attribute yield reduction to extreme weather causes, we first quantified the anomalies of precipitation and temperature.Given the rising frequency and intensity of extreme weather events, we utilized a decade's record of weather data (2013-2022) as a reference period to depict recent average conditions and identify weather anomalies.Contrasted with the commonly used 30 years period, the ten-year window better represents current occurrence of extremes allowing to more effectively inform adaptive management strategies.Marginal probability distributions have been commonly used to derive standardized anomaly indicators for meteorological variables, such as the Standardized Precipitation Index (SPI) (McKee et al., 1993) and the Standardized Temperature Index (STI) (Zscheischler et al., 2014).These anomaly indicators enable comparisons across different locations and time periods.
While SPI and STI are typically computed on a monthly basis, we used a weekly time scale for timely detection of short-term extreme weather events like heavy rainfall, heat waves, and frost.We first rescaled the gridded SILO weather data from daily to weekly for each year.Then, we fitted the weekly weather data to a Gamma distribution for each grid cell separately and subsequently standardized them by normal quantile transformation (Edwards and McKee, 1997).Since we are interested in the effects of extreme temperatures such as heat and frost, the STI was modified according to the Heat and Cold Wave Index (HCWI) to capture the extreme high and low temperatures.Specifically, we calculated the STI for both daily maximum (T max ) and minimum (T min ) temperatures, named STI max and STI min .To characterize the heat and cold waves, we replaced the threshold values from mean temperature to the 90 th percentile of T max and 10 th percentile of T min , respectively.The SPI, STI max , and STI min were calculated as: SPI = (P i − P mean )/P std (6) where P i refers to the accumulative precipitation at week i, P mean and P std refer to the average and standard deviation of precipitation for week i in the last ten years (2013-2022).T maxi and T mini represent the average of daily maximum and minimum temperature for week i.T maxQ90 , T minQ10 , T max std and T min std were calculated as the 10 years' 90 th percentile maximum temperature and 10 th percentile minimum temperature for week i and their standard deviation.For the temperature anomalies, we set as specific rules that: 1) only if the T mini goes below 0 • C (indicating frost), the calculated STI min can be negative; and 2) only if the Tmax i exceeds 37 • C (signalling a heat event that is harmful to plant growth), the STI max can be positive.
According to the SPI classification scheme used in the Copernicus European Drought Observatory (https://edo.jrc.ec.europa.eu)and the definition of STI, we defined four extreme weather indices (Eq.9-12), named the water deficit condition index (WDCI), the water surplus condition index (WSCI), the heat condition index (HCI), and the cold condition index (CCI).Increasingly severe rainfall deficits (i.e., meteorological droughts) are indicated as SPI decreases below -1.0, while increasingly severe excess rainfall is indicated as SPI increases above 1.0.The value of each index represents the intensity of potential weather extremes.

Stressors labelling
To understand which extreme events cause crop damage and to what extent, we analysed the changing velocity of the CDI curve and compared it with the four extreme weather indices (Fig. 6).We first calculated the first derivative of SG-smoothed CDI curve (CDI_1D) for each pixel.Positive derivatives indicate a period with increasing crop damage, with its peak indicating the most pronounced decline in vegetation health status.We consequently detected all peaks on the CDI_1D curve where values exceeded 0.002.The threshold was set empirically by comparing the curves and the occurrence of stressors.We also considered a second smaller peak with a CDI_1D value smaller than 0.002 but only if it occurred following the largest peak.In this way, we could capture compound effects of multiple stressors that happen intermittently during the growing season.Correspondingly, the stress events (when the extreme weather indices have a positive value) that happened within three weeks before the identified peaks were considered as the potential causes of crop damage.We applied the three-week window because it strikes a balance between sensitivity of short-term changes and relevance to crop growth dynamics, considering the variability in crop response time to extreme weather events.The three-week window provides a timeframe that is sufficiently long to detect potential events influencing crop growth variation while excluding irrelevant events that occurred one month or earlier.

Comparison of EVI2 curves from calendar and thermal time
The AGDD-based curves show a smaller deviation between the different years as compared to the DOY-based curves, especially at the early vegetative crop growth stage (Fig. 7).Wheat or barley planted in different years have similar thermal time requirements for germination and tillering.Afterwards, crops grow at different rates during the reproductive phase (stem elongation to flowering).This may be attributed to the sensitivity of crop growth to various weather extremes during the reproductive phase.

Prediction of yield reduction 4.2.1. Single variable regression
The quadratic model shows relatively higher R 2 and NRMSE between yield reduction and all individual CDI and cumCDI variables than the linear model for both wheat and barley (Fig. 8).Among the variables with strong correlation, AGDD-based variables perform better than DOY-based variables.We observed a temporal fluctuating R 2 and NRMSE among variables in different phenological stages for wheat, among which CDI in SOS90 (R 2 = 0.41, NRMSE = 0.14) and cumCDI between EOS90 and EOS70 (R 2 = 0.4, NRMSE = 0.15) show relatively higher correlations compared with other predictors.In contrast, barley shows a very high correlation during peak season, i.e., CDI in EOS90 (R 2 = 0. 69, NRMSE = 0.11) and cumCDI between PS and EOS90 (R 2 = 0.88, NRMSE = 0.08).The different correlations of wheat and barley may relate to their different responses to weather extremes.Wheat is sensitive to various extreme-weather-related stresses throughout the growing season, which can result in crop damage and yield variability occurring at different phenological stages (Cossani et al., 2009).Barley is more tolerant to stress than wheat and the causes of barley yield variability are often more concentrated and pronounced during specific growth stages, particularly during the peak growing season (Cossani et al., 2007;Dawson et al., 2015;Ryan et al., 2008).

Covariate selection
The correlation between input variables is higher when they are  S6).The most important variable for explaining wheat yield reduction is cumCDI for period EOS90-EOS70.The cumCDI between PS-EOS90 is the most important variable for predicting barley yield reduction.To reduce the potential collinearity among input variables, we marked the variables with high correlation (r ≥ 0.8, Fig. 9) as black crosses in Fig. 10a and Fig. 10b following the rule mentioned in Section 3.3.3.
Fig. 10c and Fig. 10d present the influence of the number of variables on R 2 , when adding predictors based on their average ranking result.We assessed the R 2 for the DOY-and AGDD-based time series, and for both also with and without SM transformation.The R 2 increased with the number of input variables and tended to stabilize for larger numbers of input variables.By removing the variables with high correlation (marked with black crosses in Fig. 10a-b and Fig. S6) from the selected optimal number of variables, we used the retained variables as the optimal set of predictors for further analysis (shown as green bars in Fig. 10a-b and Fig. S6).Specifically, the final selected variables include: a) for wheat, cumCDI variables between SOS50 and EOS70 for all four CDI curve construction methods, b) for barley, cumCDI between SOS50 and EOS70 derived from DOY-based methods and cumCDI between SOS90 and EOS70 derived from AGDD-based methods for barley.

Multivariate regression
We tested the yield reduction prediction ability of RF using the selected variables (Section 4.2.2) based on the DOY versus AGDD time series, with and without SMF.Fig. 11 shows the corresponding scatterplots between measured and RF-predicted yield reduction for wheat and barley over a validation set.The best performance in estimating yield reduction (R 2 = 0.83 and NRMSE = 0.09 for wheat and R 2 = 0.91 and NRMSE = 0.10 for barley) was obtained with the AGDD-EVI2 time series without applying SMF.For both wheat and barley, AGDD substantially improves the prediction accuracy compared with the DOYbased method, whereas the combination of AGDD and SMF method does not perform better than solely using AGDD.The SMF method applied for DOY curves improves the prediction performance, but not as much as AGDD transformation.When comparing the distribution of RF-predicted yield reduction against measured yield reduction (Fig. 12), the model tended to overpredict small yield reduction and underpredict large yield reduction using all four CDI construction methods, especially for wheat.This discrepancy may stem from a combination of factors, such as differences in crop sensitivity and resilience to stress between wheat and barley, farmers' remedial interventions for unhealthy crops, and uncertainties in reference curve construction and CDI calculation.
The comparison results of different VIs, CDI curve construction methods, and regression models are shown in Fig. S7 and Fig. S8.Among the four VIs, EVI2 performed better than the other three VIs tested.The RF model consistently outperformed support vector regression and multiple linear regression, yielding higher R 2 values in most instances.The prediction accuracy differed substantially between DOY-and AGDD-based time series, depending on the regression model utilized.Overall, the AGDD-EVI2 time series, without SMF applied, exhibits the highest accuracy in estimating yield loss for both winter wheat and barley when employing the RF model.

Timing and attribution of weather extremes
Fig. 13 shows an example for a single pixel of the interannual variations in detected extreme weather events that may be linked to the timing of crop damage.The attribution analysis of crop damage is done separately for each pixel within the studied paddocks following the stressors labelling rules explained in Section 3.4.2.With the detected peaks (> 0.02) in CDI_1D curve (i.e., the first derivative of CDI curve), a three-week window before this peak is applied to define a period of poor crop development indicating the likely occurrence of crop yield reduction.This period is then compared with the derived time series of the extreme weather indices (HCI, CCI, WDCI, and WSCI) to evaluate if the poor crop development coincides with a weather event that occurred during that period.
Fig. 14 presents the spatial and temporal variation of within-paddock crop yield reduction and its potential attribution for three years taking one paddock as an example.This paddock was cultivated with wheat in 2017, 2020, and 2021.The yield reduction was estimated at the pixel level by applying the RF regressor based on cumCDI variables between SOS50 and EOS70 extracted from the AGDD-EVI2 time series (i.e., the optimal phenoperiod for yield reduction estimation selected in Section 4.2).Crop damage in 2017 was mainly caused by water deficit and frost that occurred in July.The stressor map for 2020 shows that the spatial variation of crop damage may be attributed to different weather extremes even within the same paddock, including frost in August and drought in July; the strongest yield reduction was found for areas that experienced a combination of frost and drought.Wheat yield reduction is smaller in 2021 as compared with other years.Water surplus in August is the dominant reason of crop yield reduction in 2021.The observation of multiple stressors within a single field may be explained by the following reasons: 1) for extreme precipitation small-scale variability can be substantial even if not captured by our rainfall data, 2) soil properties and micro-topography differ within a field leading to spatial variability in how much water is received or retained.
Finally, we produced the crop yield reduction maps and corresponding attribution maps of weather extremes for all wheat and barley paddocks in 2017-2022 (Fig. 15 and Fig. 16).The results were validated by comparing against auxiliary data (Table S2) and showed consistency with reports on the actual stress.For example, stem frost impacted crop yield in 2020 in most wheat cultivation areas while only a few barley paddocks suffered from frost damage.Wheat is more susceptible to frost than barley at the flowering stage (Ferrante et al., 2021).Due to a rainfall gradient in the study area that increases from north to south, the northern part usually experiences more intense drought conditions.This pattern is confirmed by the yield reduction maps for barley in 2019 and 2021 which show a higher yield reduction in the north.Fig. 8. R 2 and NRMSE of each CDI and cumCDI variable for predicting wheat and barley yield reduction using linear and quadratic regression based on DOY-and AGDD-EVI2 time series.The labelling around the radar chart refers to the abbreviated CDI and cumCDI variables (e.g., "S10″ means "CDI_SOS10", "cS10-S30" means "cumCDI_SOS10-SOS30").

Advantages of the proposed approach
We proposed a scalable system for early prediction of winter cereals yield reduction within the growing season and for attributing this reduction to the occurrence of extreme weather events at sub-field level.The method incorporates a CDI using EVI2 time series derived from Sentinel-2 data that are scaled to AGDD based on weather information.By analysing the CDI curve within specific phenological windows, we demonstrated that CDI information between the stem elongation and grain filling stage can accurately predict potential yield reduction.Furthermore, by comparing moments of significant CDI increase with different indices on weather extremes that coincide with those moments, we identified the potential extreme weather events responsible for the observed yield reduction and determined the timing of these stresses.
In contrast to directly using EVI2 time series, our proposed CDI transformation shows a significant improvement (R 2 increased by 0.31 for wheat and 0.13 for barley compared to the baseline) in yield loss estimation.We conducted a comparative analysis with a baseline model that employs the RF algorithms with DOY-and AGDD-based EVI2 time series as direct inputs.Instead of detecting phenometrics, we derived time series of monthly average EVI2 as DOY-based features and applied EVI2 values at steps of 300 degree days as AGDD-based inputs for the RF model.The results (Fig. S9) of this comparison underscore the superiority of CDI features, especially AGDD-based CDI variables without SMF (with high R 2 values of 0.83 for wheat and 0.91 for barley), in predicting yield loss.The predictions of the baseline model exhibit the lowest R 2 (0.52 for wheat and 0.78 for barley) and significant differences (p ≤ 0.05) compared with both CDI and SMF-fitted CDI models in wheat and barley.The discrepancy highlights the effectiveness of our proposed CDI phenometrics-based model in reflecting the physiological status of the crops under stress and representing the complex relationships between the influence of stress and crop damage.
We evaluated the performance of different profile construction methods and phenological phases for the prediction of yield losses at the sub-field level.Our results agree with results from previous yield estimation studies while also suggesting important improvements.Firstly, Hunt et al., 2019 demonstrated that RF performed better for within-field yield estimation as compared to more traditional regression models (e. g., linear regression model), which is in line with our findings.The integration of detected sowing dates to recalibrate the start point of the VI time series was also applied by Funk and Budde, 2009 who used this to improve yield prediction accuracy.Our results enhanced the yield loss estimation accuracy by transforming DOY to AGDD-based time series using a similar phenology-adjustment method as Funk and Budde, 2009.Furthermore, previous research indicated that the accurate yield prediction can be achieved two-month before maturity and the optimum phenological window concentrated on the mid-to-late season (Benedetti and Rossini, 1993;Cai et al., 2019;Doraiswamy and Cook, 1995;Funk and Budde, 2009).Similarly, our results underscore the superior predictive capability of the period between SOS50 and EOS70, corresponding to the stem elongation and grain-filling stage.
In contrast to existing research focused on yield loss attribution analysis at post-season based on crop model and statistical algorithms, our approach applied high-resolution VI time series, offering the advantage of timely identification of crop damage with respect to the occurrence and timing of weather extremes.Prior efforts calculated the contribution of diverse extreme-weather-related stresses to yield reduction utilizing regression models or crop simulation models (Nóia, 2023;Vogel et al., 2019;Zhao et al., 2022;Zhong et al., 2023).However, not all instances of extreme weather result in yield reduction, depending on both the timing of occurrence and crop susceptibility.Furthermore, such analysis is typically conducted post-season or late-season due to the required yield data. Solely relying on the correlation between extreme weather indices and yield data might prove insufficient in accurately  attributing yield reduction to its causal factors.We addressed these gaps by introducing an innovative method that elucidates the within-field variability in yield reduction and its attribution, facilitating the identification and comprehension of the impacts of various extreme weather events on crop health at the sub-field level.
Incorporating thermal time in the EVI2 time series proved beneficial to diminish temporal deviation between time series across multiple years, contributing to the construction of typical stress-free crop growth curves, better detection of crop anomalies, and prediction of yield losses.Nyborg, Pelletier, and Assent (2022) reported that using thermal-timebased time series improved crop classification by considering varying crop growth rates under different years and weather patterns and thus reducing the temporal shift.Moreover, Zeng et al. (2016) demonstrated that the use of thermal time rather than calendar days more accurately  Wheat yield reduction maps and corresponding spatial distribution and timing of attribution to weather extremes from 2017 to 2022 for all paddocks that were available in our dataset.Grey polygons correspond to paddocks that were not cultivated with wheat or without crop type information in the respective year.Letters F, D, and W in the corresponding legend represent frost, drought (water deficit), and water surplus, respectively.captures phenology variability and may enhance yield prediction, considering the close relationship between crop yield and critical phenology stages.We confirmed this assertion, as the prediction of yield reduction applying thermal-time-based (i.e., AGDD-based) EVI2 time series outperformed calendar-based curves.This highlights the importance of factoring in temperature, a primary environmental factor controlling the growth of winter cereals, when developing crop loss detection models (Slafer and Rawson, 1994).

Scope for improvement
Despite the good performance of thermal time in diminishing differences in crop growth rates at early growing stages, other environmental factors such as water availability can also have large effects on crop development.Depending on climate and crop, cropland productivity can be temperature-limited (Akter and Rafiqul Islam, 2017;Xiao et al., 2018), moisture-limited (Blum, 1996;Brown and de Beurs, 2008;Li et al., 2009), radiation-limited (Nemani et al., 2003), or affected simultaneously by multiple limitations (Pradhan et al., 2012;Shah and Paulsen, 2003).Identifying and integrating the appropriate environmental variables, or a combination thereof, is critical to accurately model crop growth status and compare this status more effectively between locations and years.For example, accumulated relative humidity can be used as an alternative to thermal-time in RS time series for semiarid areas where crop growth is mainly moisture-limited (Brown and de Beurs, 2008).Zeng et al., 2016 proposed the concept of photothermal time combining both temperature and photoperiod information and illustrated the potential of minimizing inter-annual variation of detected phenological stages in irrigated areas.This concept could be transferred and assessed within our method.Nevertheless, challenges still remain in defining a variable that quantifies the combined relationship between temperature, water, and radiation, which can be incorporated into the RS time series to improve model transferability across different crops and climate conditions.
Another limitation is that the used rainfall and temperature data may not accurately represent fine-scale differences between and within paddocks in our study area.The gridded weather data with five km resolution can overlook micro-climate variations, which are usually influenced by local topography, soil properties, vegetation cover, and agricultural practices.As a consequence, the extreme weather indices that we derived may not accurately represent the exact stress events experienced at specific locations, possibly leading to an erroneous attribution of yield reduction if using such weather-based stress indicators.To illustrate, the southern regions of our study area feature undulating topography, contrasting with the relatively flat terrain in the north.This divergence influences rainfall distribution and contributes to more pronounced drought conditions in the north.While the yield reduction maps capture such spatial gradient, the stress maps fail to mirror this variation due to the spatial uniformity of SILO data across the study region.Recognizing this limitation underscores the necessity for more precise and detailed sources of meteorological data capable of accounting for small-scale weather variability as caused for example by elevation.For instance, supplementing the existing network of weather stations with terrestrial radar data or leveraging microwave links between antennas of mobile networks could yield more accurate and localized insights on precipitation distribution (Graf et al., 2020;Overeem et al., 2013).
In addition to the attribution of crop yield reduction to extreme weather events, it is important to acknowledge the potential influence of other biotic and abiotic factors on crop health and productivity.While our study focused on identifying whether extreme temperature and precipitation events could partially explain yield loss, crop losses can result from a combination of factors beyond the scope of our analysis.Pest infestation, diseases, soil quality, management practices, and socioeconomic factors all play crucial roles in shaping crop health and productivity (Gull et al., 2019).Moreover, these factors often interact in complex ways, amplifying the challenges of accurately attributing yield loss to specific causes.While extreme weather events can have significant impacts on crop yields, their effects may be compounded or moderated by other stressors (Suzuki et al., 2014).Therefore, the mere attribution to extreme weather events is a simplification of the complex reality of agricultural systems.Moving forward, a more holistic approach that considers the multifactorial nature of crop yield reduction could be considered by integrating information on various stressors and their interactions.However, such information may be difficult to obtain across wide geographies.While a more comprehensive understanding of crop loss dynamics could inform more effective mitigation and adaptation strategies in the face of changing environmental conditions, understanding to what extent extreme weather events link to crop loss is a first important step to understand and adapt to climate challenges in agriculture.

Pathways for future development
While initially tested on winter cereals in South Australia, our method holds potential for extension to other regions and crops.Future research should prioritize enhancing the model's adaptability through testing it with a wider range of crop datasets and exploring transfer learning strategies to improve its transferability in space and time.However, challenges and limitations may arise due to the diversity of crop types, regional agricultural practices, and meteorological conditions.
In addition to the model improvements detailed in Section 5.2, which involve the exploration of additional weather-based variables and more detailed meteorological data sources, careful attention should be given to constructing the reference curve and addressing potential data quality issues.Expanding the application of CDI curves to larger regions necessitates redefining the spatial scope for establishing reference curves.Single reference curves were developed for both wheat and barley considering the small spatial extent of the study area.However, when expanding the analysis to larger regions (e.g., the country or continental scale), the utilization of a single curve to represent the typical growth trajectory of a crop becomes inappropriate due to the diverse environmental conditions and distinct crop varieties.A solution could be to calculate a reference curve for each small region (e.g., at the county level) based on crop type maps in the historical years or alternatively using relatively homogenous zones as determined by RS-based ecological stratification (De Oto et al., 2019).Moreover, when applying CDI curves to new areas, data gaps and quality issues with optical image time series may arise.While we used Sentinel-2 data, which offers global coverage with a revisit time of five days in principle, cloud cover can lead to missing data and low-quality observations.Strategies such as data fusion methods, which involve combining different sources of satellite images to fill in missing data gaps, can help mitigate this limitation.For example, Synthetic Aperture Radar (SAR) data, like Sentinel-1, offer high resolution images immune to cloud cover or atmospheric interference.These SAR images can complement Sentinel-2 data, filling gaps and ensuring continuous coverage of the study area (Garioud et al., 2021;Pipia et al., 2019).Another example is the Harmonized Landsat Sentinel-2 (HLS) product, which combines observations from both Landsat and Sentinel-2 satellites, enhancing temporal coverage and consistency (Claverie et al., 2018).Addressing these aspects is crucial to ensure the scalability and effectiveness of our approach in detecting and attributing crop yield losses accurately and reliably across diverse agricultural landscapes.
Our framework (i.e., crop damage detection and attribution with CDI) provides a basis for improving the future prediction capability of crop yield reduction with RS time series, for early warning of crop damage, and for attributing its weather-related causes in near real-time, which can help, among others, to adopt appropriate crop protection strategies.With the assistance of high-resolution Sentinel-2 data and daily meteorological data, the CDI curve and stress indices can be calculated in near real-time.Predicting potential yield losses using AGDD-based time series involves identifying the onset date, which using the methods of this paper can be derived once the peak season (i.e., maximum EVI2) occurs.Subsequently, yield reduction can be estimated based on the cumCDI between SOS50 and the latest-available observation given that SOS50 to EOS70 is the optimal period for yield loss prediction.Correspondingly, potential yield reduction can be predicted at least two months prior to harvesting, specifically during the stem elongation to grain filling stage.Crop damage can be detected once a small peak on the first derivative of CDI curve is observed; the damage can then be attributed to specific extreme weather events as CDI increases peak usually within three weeks after the extreme events occur.However, the timeliness may be compromised by the need for temporal filtering.For example, the SG filter smooths raw time series and requires data available before and after each observation with a specified window size, aiming to remove noise and capture underlying trends.The choice of window size can impact the trade-off between time series smoothing and data details preservation, thereby affecting the timeliness of crop damage detection.Different from the commonly used smoothing method, near real-time filtering provides higher flexibility by projecting to the current time only from past observations (Sedano et al., 2014).The reliability of filtered estimation at each time step improved with additional observations becoming available.However, the estimation errors especially at the initial consolidation stage may be amplified by further anomaly detection (Meroni et al., 2019).To improve product quality, the use of satellite data with high revisit intervals (e.g., PlanetScope) could be a solution to improve both the accuracy and timeliness for near real-time filtering.
In summary, we conducted a proof-of-concept study and demonstrated the potential of the proposed framework for near real-time monitoring of crop yield losses and their accurate attribution to weather extremes.This initial investigation sets the stage for future studies, which may include detailed assessments aligning with microclimate conditions or more comprehensive evaluations covering larger geographical areas.As such, it can play a crucial role in the framework of climate-smart agriculture given that yield loss and attribution information aids in proactive adaptation to a changing climate by enabling informed decision-making and responsive strategies for sustainable agricultural practices.

Conclusion
We presented a novel and scalable method for predicting in-season yield reduction at the sub-field level and attributing these losses to extreme weather events.Using temporal profiles of a crop damage index (CDI) derived from thermal-time-based VI time series, we could accurately predict yield reduction of winter cereals, and attribute this reduction to the likely extreme weather causes.The use of thermal-timebased VI time series reduced temporal deviation of EVI2 curves of different years and locations, and by better aligning the crop growth development, it outperformed calendar-based curves in yield reduction prediction.Yield reduction could be accurately predicted for wheat and barley approximately two months before harvest, with high R 2 values of 0.83 and 0.91, respectively.The combined analysis of CDI curves and extreme weather indicators enabled the timely detection of weatherrelated causes of crop damage at the sub-field level, providing valuable insights for improving crop management responses.
This study stands as a proof-of-concept by testing the method for winter cereals in South Australia.It provides a basis for early warning of crop damage and attributing such damage to weather extremes in near

Fig. 1 .
Fig. 1.Overview of the study area and crop paddocks (black polygons) for which cropland information was provided by local growers for 2017-2022, with an example of winter wheat and barley paddocks in 2021.Fig. S1 in the supplementary materials shows the crop maps for the other years.

Fig. 2 .
Fig. 2. Growth stages of winter cereals (wheat and barley) in Australia based on the Zadoks decimal growth scale.Adapted from J. Hunt et al., 2018.

Fig. 3 .
Fig. 3. Weekly (a) maximum temperature (maxT), minimum temperature (minT), and (b) precipitation (P) acquired from the SILO weather database for 2017-2022.The black curves represent the average value (ave) for the last 10 years (2013-2022).The colours refer to data corresponding to different years from 2017 to 2022.

Fig. 4 .
Fig. 4. Flowchart of the analysis.Numbers in dark blue (n=…) refer to number of in-situ samples (or directly derived from these samples) for wheat and barley paddock years combined.Abbreviations: EVI2 − two-band enhanced vegetation index, DOY − day of year, AGDD − accumulated growing degree days, CDI − crop damage index.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Representation of the CDI curve derivation.The steps include: (a) establish the reference curve based on historical data, (b) extract phenometrics from the observed curve, (c) fit the reference curve (i.e., shape model) to the observed curve at vegetative stage, and (d) construct CDI curve by calculating the differences between scaled shape model and observed curve.Abbreviations: SOS − start of season, EOS − end of season, PS − peak season.The numbers behind refer to the percentage of the difference between the maximum and minimum EVI2 in upward and downward phases.

Fig. 6 .
Fig. 6.Schematic illustrating the fundamental principle employed for detecting weather-related factors attributing to crop damage.The blue bars indicate the detected crop damage causes (i.e., extreme weather events which happened within three weeks before the identified peaks in CDI_1D curve).Abbreviations: CDI_1D − first derivative of CDI curve, HCI − heat condition index, CCI − cold condition index, WDCI − water deficit condition index, WSCI − water surplus condition index.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7.The DOY-and AGDD-based EVI2 time series of winter wheat and barley paddocks for 2017 to 2022.Each colour represents the average time series of all pixels for the available paddocks of target crop in one year.The circle represents the observation data and the curve refers to the SG-filtered time series.

Fig. 9 .
Fig. 9. Pearson correlation coefficient matrix (significant at p ≤ 0.01) between variable pairs derived from wheat (a) DOY-EVI2 and (b) AGDD-EVI2 time series and barley (c) DOY-EVI2 and (d) AGDD-EVI2 time series.Positive correlations are displayed in red and negative correlations in blue.Colour intensity means proportional to the correlation coefficients.White spaces refer to the variable pairs for which the correlation was insignificant (p > 0.01).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.(a-b) Feature importance value rank for wheat and barley based on variables derived from AGDD-EVI2 and (c-d) the corresponding change of R 2 (based on testing dataset) for the prediction of yield reduction with an increasing number of input variables gradually based on the feature importance rank.The black cross behind each bar in panels (a-b) represents the variable which has a high correlation (r ≥ 0.8) with another variable at earlier phenological stage.The features for which a green bar is shown are the final retained variables for multivariate regression analysis.The vertical dash lines in panel (c-d) represent the selected number of variables where the curve stabilizes (i.e., R 2 increase from n to n+1 is less than 5%).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Raincloud plots illustrate the distribution of in-situ yield reduction records and predicted yield reduction using different CDI construction methods (DOY, AGDD, DOY-SMF, AGDD-SMF) based on the RF model for wheat and barley.The height of each violin represents the density of yield reduction at different values along the x-axis.Under each violin, scatter plots provide a detailed view of individual data points.The box plot denotes the median value of yield reduction and its interquartile range (25-75%), and the whiskers extending from the box refer to the rest of the distribution.

Fig. 13 .
Fig. 13.Illustration of the detection of the timing of wheat yield reduction and its attribution to weather extremes.The CDI curves show an example pixel of wheat paddock.The orange shade in the two top rows indicates the crop growing season.The stress events within the growing season are highlighted in orange in the four bottom rows.The red line in the second row refers to threshold of 0.002 for detecting peaks in CDI_1D curve.The green box represents the 3-week windows around the detected peaks.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .
Fig. 14. (a-c) Yield reduction map derived from RF model applied to the AGDD-EVI2 time series and (d-f) potential stressor map of an example paddock cultivated by wheat in 2017, 2020, and 2021.The different colours in (d-f) refer to different kinds of weather extremes.Letters F, D, and W in the corresponding legend represent frost, drought (water deficit), and water surplus, respectively.

Fig
Fig.15.Wheat yield reduction maps and corresponding spatial distribution and timing of attribution to weather extremes from 2017 to 2022 for all paddocks that were available in our dataset.Grey polygons correspond to paddocks that were not cultivated with wheat or without crop type information in the respective year.Letters F, D, and W in the corresponding legend represent frost, drought (water deficit), and water surplus, respectively.

Fig. 16 .
Fig. 16.Barley yield reduction maps and corresponding spatial distribution and timing of attribution to weather extremes from 2017 to 2022.Grey polygons correspond to paddocks that were not cultivated with barley or without crop type information in the respective year.Letters F, D, and W in the corresponding legend represent frost, drought (water deficit), and water surplus, respectively.