The region of interest is situated in north-eastern Germany within the federal state of Mecklenburg West Pomerania (Figure 1). Characterized by an average annual precipitation of 550 mm and a mean air temperature of 8.3 °C, its climate is classified as temperate Middle-European. By 2001 the German Aerospace Center (DLR) established the Durable Environmental Multidisciplinary Monitoring Net Information Network (DEMMIN), which serves as a calibration and validation site for Earth observation missions (Spengler et al., 2018). DEMMIN is also part of the Joint Experiment for Crop Assessment and Monitoring (JECAM) network (Hosseini et al., 2021).

For the growing seasons from 2017 to 2021, information on parcel delineation and crop types in the DEMMIN area was extracted from the German integrated administration and control system (InVeKoS) of Mecklenburg West Pomerania. Winter wheat, canola and sugar beet are major crop types in DEMMIN. They belong to the seven area dominating crop types in Germany and were hence selected for this study. Only fields of significant size (>3 ha) were chosen to minimize pixel contamination from neighbouring land use or cover as previously described (Lobert et al., 2023; Löw et al., 2024).

Phenological in situ observations were obtained from the voluntary-based observer framework of the DWD (Kaspar et al., 2015). Any phenological developments recorded by DWD were translated to the BBCH (Biologische Bundesanstalt für Land-und Forstwirtschaft, Bundessortenamt und CHemische Industrie) scale, i.e., a classification scheme for measuring phenological development (Meier, 2001) to facilitate comparison with other studies. These observations encompass emergence, leaf development and canopy closure of sugar beet, stem elongation (BBCH 30), heading (BBCH 50), yellow ripening (BBCH 87) and harvest (BBCH 99) for winter wheat as well as inflorescence (BBCH 50), start (BBCH 60) and end of flowering (BBCH 69) and harvest (BBCH 99) of canola. This dataset was used to analyse approximately 500 field of winter wheat, 300 fields of canola and 150 fields of sugar beet per year. Due to the area’s poor coverage by this monitoring setup in terms of proximity and crop variety, the average occurrence date by federal state served as source of validation.

Growing degrees are characterized as heat units and serve as a widely utilized tool for elucidating the progression of biological processes (McMaster and Wilhelm, 1997). Growing degree days (GDD) represent the cumulative sum of these units employed to model crop growth. They were also used in remote-sensing based studies to translate crop phenology into crop maturity for integration into crop yield models (McNairn et al., 2018). However, classical GDD approaches rely on a minimum base temperature only and neglect an upper temperature limit for plant growth. Especially conditions of extreme temperatures and variability can lead to errors when using those GDD approaches for phenological assessments (Ritchie and Nesmith, 2015; Stinner et al., 1974; Zhou and Wang, 2018).

Given our observation period, which includes years marked by extreme drought in Germany (Harfenmeister et al., 2021b; Schlund and Erasmi, 2020) the method proposed by Zhou and Wang (2018) was selected for meteorological assessments of GDD. This method addresses the shortcoming by the calculation of hourly temperature time (HTT) series with four key parameters: T _h (hourly measurement of air temperature), T _u (upper temperature limit), Topt (optimal temperature for maximum growth), and T _b (base temperature) (Zhou and Wang, 2018). By implementing an optimal threshold and an upper temperature limit to HTT calculation it can reflect heat and other temperature related stress, because less or no additional HTT are accumulated (see Equation 1). H T T = 0 , T h < T b T h − T b T o p t − T b T u − T h T u − T o p t T u − T o p t T o p t − T b   T b ≤ T h ≤ T u 0 , T u < T h (1)

Data on air temperature is provided by DWD for the entire observation period (2017–2021) as an interpolated raster derived from the sensor network of DEMMIN (Haßelbusch and Lucas-Mofat, 2021). The native spatial resolution of the dataset is 250 m × 250 m per pixel, its temporal resolution is one hour and it was interpolated by employing ordinary kriging. Further details can be found in this DWD-report (Haßelbusch and Lucas-Mofat, 2021). Thus, the GDD are then calculated as cumulative sum of daily mean HTT based on this dataset. These calculations were conducted on an open data cube platform (ODC) (Killough, 2018) using the programming language Python. Table 1 specifies the thresholds of air temperature by crop type according to named references.

We used values derived from literature, because no data on strains and phenotypes for the investigated crop was available for that time period and area.

The time series of S1 data encompasses relative orbits 146, 95, and 168, spanning the period from 2017 to 2021, amounting to approximately 1.050 data sets. S1 data sets were obtained in Interferometric Wide Swath mode and at VV/VH polarization. Table 2 lists the pass direction as well as the range of incidence angles for each relative orbit.

The calculation of interferometric and polarimetric features necessitated the use of Single Look Complex (SLC) data (ESA, 2013). The data pre-processing involved a combined approach utilizing pyroSAR (Truckenbrodt et al., 2019), a Python-based API, and SNAP (Version 9) (ESA, 2025). This processing framework was seamlessly integrated into an ODC environment, which served as the data management platform. This platform runs on a cluster located at the Leibniz-DataCentre (LRZ) and is equipped with sufficient computation power to process at tile of SLC data in a few minutes (Friedrich et al., 2024). The polarimetric feature processing sequence comprised terrain flattening (Small, 2011), multi-looking with one look in azimuth and four looks in range, speckle filtering through a 5 pixel x 5 pixel boxcar filter, and Range-Doppler Terrain correction (Richards, 2009), yielding a spatial resolution of 20 m × 20 m and gamma nought (GN) backscatter. Given the dB scaling of backscatter, the cross-pol ratio (CR) was computed as VV-VH. Alpha and Entropy were derived from a C-2 Matrix (Cloude and Pottier, 1996). This set of S1 features consisting of VV/VH backscatter, CR, Alpha and Entropy as well as VV/VH coherence, has been proven to reflect the changes in plant physiognomy along the crop life cycle by various studies. No further decomposition metrics were included, because they possess less explanatory power in the dual-pol data of S1 when compared to quad-pol data from e.g., Radarsat. (Harfenmeister et al., 2021b; Khabbazan et al., 2019; Löw et al., 2021; Meroni et al., 2021; Schlund and Erasmi, 2020).

For interferometric (InSAR) coherence, the parameters for multi-looking and terrain correction remained consistent. Coherence calculation involved removing the flat earth and topographical phase using a moving window of three pixels in azimuth and eleven pixels in range, along with a six-day temporal baseline and consecutive reference images (Zhang et al., 2018). SRTM data (1 Arc-second) served as the digital elevation model in all steps requiring it. This set of S1 features, comprising VV/VH backscatter, CR, Alpha, and Entropy, as well as VV/VH coherence, demonstrated its ability (Harfenmeister et al., 2021a; Khabbazan et al., 2019; Löw et al., 2024; Schlund and Erasmi, 2020) to effectively capture changes in plant physiognomy related to phenological developments in various studies.

Since this study also investigated potential explanations of the spatiotemporal distribution of average agreement (AVA), dominance of tendencies (DoT) and outlier occurrences, auxiliary data was acquired and processed. The collected data represents various environmental factors that may explain the occurrence of these metrics. In regard to soil properties SoilGrid (Poggio et al., 2021) data on soil organic content (SOC), organic carbon stock (OCS) and the dominant soil type was downloaded for the area. SOC and OCS are used as proxy for soil fertility in this study. It is assumed that plants display a different physiognomy depending on the fertility of the soil. Because SAR based phenology is sensitive to changes in plant physiognomy SOC and OCS have an indirect impact on the interaction of signal and plant physiognomy. Furthermore, SOC and OCS are related to soil texture and therefore surface roughness, which is also an influencing factor in early phenological stages, when there is still soil signal present (Santos et al., 2023). In addition, daily sums of precipitation and mean air temperature were calculated from the same interpolated dataset that was used to calculate GDD (Haßelbusch and Lucas-Mofat, 2021) to support the discussion of the results in terms of local water availability. Moreover, terrain features such as aspect, slope as well as the topographical wetness index (TWI) (Sørensen et al., 2006) were derived from an SRTM30 digital elevation model. These features plus the elevation were used to explain distributions by the terrain itself. Lastly, pixel-based, mean standard deviation of yearly time series of enhanced vegetation Index (EVI) and normalized difference vegetation index (NDVI) were used as explanatory variables representing vegetation dynamics and field heterogeneity (Gessner et al., 2023). Both Indices were calculated using Sentinel-2 L2A data that were cloud-masked by the scene classification layer (SCL).

The idea behind the following approach is a comparison of various time windows. For example: In 2020, at landscape level, the TSM distribution of VH intensity of canola level displayed many relevant TSM occurrences around day of year (DOY) 115. By analysing single fields it can now be ascertained which canola fields show a similar behaviour to that pattern observed at the landscape scale. Due to the GDD baseline and distances between associated GDD values the similarity or the lack of it can be quantified. Similarly, these findings were linked to the in situ observations made by DWD. This process is applied to each S1 feature, orbit and for every year that is covered by the catalogue of Löw et al. (2024).

All indicators can be calculated for each year, relative orbit, S1 feature and crop type. Moreover, field specific uncertainties can be computed for each BBCH stage. Depending on the required level of detail, these indicators can be aggregated accordingly. For this study, the highest level of detail is year, crop and S1 feature specific, as the impact of relative orbits was already extensively analysed by Löw et al. (2024). Because the study at hand and Löw et al. (2024) share the same database a repetition of analyses identifying suitable combinations of relative orbits and S1 features was deemed unnecessary.

Subsequently, a conceptual approach was developed that assesses the spatial distribution of the above listed indicators. This assessment addresses the explainability of said indicators via the geographical settings of the study area. In an exemplary analysis, the plausibility of outliers at field level were investigated via Spearman’s rank correlation coefficient and the variable importance of a random forest (Breiman, 2001; Conrad et al., 2017). It was hypothesized that deviations at the field level from the landscape pattern depend on the field-specific environmental conditions. Hereby, a variety of environmental descriptors were included. An exemplary set of indicators is illustrated in Figure 3. The full list and the corresponding abbreviations are listed in the Supplementary Appendix (see Supplementary Table SA1). Matching pairs for the correlation analysis and training the random forest were extracted as mean values per field, except for soil type, where the dominant class was derived by modus.

The internal model accuracy of the various random forests was evaluated using root mean square error (RMSE, mean absolute error (MAE) and R squared. These metrics were derived by the internal accuracy assessment provided by the R package caret (Kuhn, 2008).

Figure 4 presents an example of field specific information on winter wheat across the years 2017, 2018, 2020 and 2021, illustrating the mean agreement across all relevant S1 features. The substantial proportion of fields exhibiting an AVA below 40% in 2017 and 2021 reflects the findings of the yearly, S1 feature specific analysis (see below). However, the majority of fields across all years were allocated within the 50%–80% agreement range. Also, fields of higher rates of agreement were detected in the south eastern, south western and north eastern region across all depicted years. Notably, fields exceeding 90% of AVA concentrated in 2020, suggesting a greater homogeneity in crop development. According to yearly, S1 feature specific statistics, winter wheat exhibited comparatively low rates of agreement, particularly for VV and VH intensity in the years of 2017 (40%) and 2021 (60%). While canola displayed a similar gradient when comparing VV and VH intensities with other S1 features, the yearly performance is lowest in 2018 and 2019 Nevertheless, canola’s overall the agreement rate is greater than winter wheat, but lower than sugar beet. In contrast, sugar beet consistently achieved the highest agreement rates, exceeding 85%, whereas canola and winter wheat mainly ranged around 70%.

As shown in Equation 3, the information on AVA was initially derived for each year and S1 feature. This provided spatially explicit information on the extent to which each field aligns with the TSM and GDD occurrences at landscape level.

The trackable progress was calculated for each BBCH stage, relevant S1 feature and year within the observation period. This resulted in large quantities of vector data, exceeding the scope of this paper for comprehensive presentation. Therefore, Figure 5 only displays an exemplary data set derived from VV coherence for wheat fields in 2017.

The left panel of Figure 5 illustrates the minimum tracking progress of GDD coinciding with the first relevant TSM occurrence in the S1 signal. The large number of fields fall within the range of 85%–90%, which is congruent with the DoT findings, where most winter wheat fields exhibit advanced development relative to the time window at landscape level. However, a subset of fields falls within the range of 110%–115%, indicating later development. The right panel of Figure 5 displays the range of GDD-based progression covered by TSM occurrences. In this example, the fewest fields were observed within the range of 30%–34.5%.

By plotting the field and BBCH specific tracking range for the same S1 feature across multiple years (see Figure 6), it became evident that there is a notable variance in the distribution of values. The years 2017 and 2021 depict a higher number of fields within larger ranges (20%–35%) whereas 2018 and especially 2019 were dominated by fields with lower ranges (0%–10%). Furthermore, distinct spatial cluster of consistently high or low ranges could be identified.

By inspecting an exemplary field specific timeline (see Figure 7), a semi-regular pattern emerged. Relevant occurrences of break points were tracked for the selected field around Julian Days 50, 100 to 125 and 175 to 225. Additionally, the distribution of TSM occurrences appeared more compact in years with higher rainfall intensities, which aligns with a larger number of fields exhibiting wider tracking ranges in 2017 and 2021 (see Figure 6). Furthermore, rainfall events did not regularly overlay with the occurrences of break points. Besides, when comparing the EVI time series with the occurrences of break points, the observed patterns are mostly congruent with the overall trends depicted by the EVI time series, despite some inconsistencies in the data availability.

At field level, detailed information about DoT and its spatial distribution is displayed (see Figure 8). The distribution of time windows appears skewed, which is in line with the findings observed at the landscape scale, which are described in a subsequent paragraph. Furthermore, information on the strength of the DoT was illustrated. In 2017 (see Figure 8), for example, sugar beet predominantly depicts fields of strong positive consistencies, whereas in 2020 and 2021 a higher proportion of fields displayed negative DoT, consistent with the yearly, S1-specific observations. Similarly, the occurrences of 2018 are consistent with the dominance of negative rates in three S1 features. In this map no clear spatial clustering of DoT across multiple years was found. This contrasts the findings in the example of AVA (see Figure 4), where clusters of similar agreement rates could be found across multiple years.

The yearly, S1 feature specific analysis of DoT revealed, that most fields of winter wheat and canola (see Figure A6) display a dominant positive DoT across all S1 features and years. In contrast, sugar beet depicted a more variable pattern. For example, coherences in 2017 were mostly dominated by the equal category. In 2018, the behind category was dominant for Alpha, Entropy and CR, whereas similar distributions of categories were observed only for VH coherence and VV intensity in 2019.

The analysis of deviations and outliers revealed that the individual field deviated from the landscape pattern up to 200 times (see Figure 9). However, these numbers reflect the total count across all observed years, crop types, S1 features, and relative orbits. By examining Figure 9, it becomes evident, that the majority of fields displays between 20 and 80 outliers. Higher outlier counts (above 120) are mostly concentrated within a few spatial clusters, particularly in the North-Northeast and the Southwestern part of the study area.

As addressed in chapter 2.7, an analysis of potential environmental explanations was conducted. The correlation analysis yielded weak, but significant (p = 0.05) correlation coefficients for the following variables elevation (elev_mean; 0.26), dominant soil type (soil_type 0.26), soil organic content at 0–5 cm (soc_mean; −0.23) and organic carbon stock at 0–30 cm (ocs_mean; −0.23). These results are visualised in Figure 10. The strongest correlation was found between organic carbon stock (ocs_mean) and elevation (elev_mean) (−0.45).

In accordance with the correlation analysis, the variable importance of the random forest revealed elevation as the most important predictor of the total outlier count. The second and third most important variables were organic carbon stock at 0–30 cm (ocs_mean) and soil organic carbon content at 0–5 cm (soc_mean) (Figure 11).

The accuracy assessment yielded the following results: RMSE = 33.5, MAE = 26.2 and R² = 0.25. Additionally, the crop specific analysis of outlier occurrences and deviations from the landscape pattern provided the following insights, summarised in Table 4.

The model accuracy for all three crop types was consistent with values around 0.85. Soil organic carbon content (SOC) exhibited a weak correlation with the outlier occurrences of all crop types. Furthermore, terrain features such as elevation and TWI were ranked as the variables of highest importance.

The exemplary analysis of the AVA showed that fields tend to agree with the overall pattern at landscape level between 60% and 70%. Hence it is assumed that the predefined threshold of 70% is valid for designating fields or in a wider context also years and S1 features of high agreement. Moreover, the indicator demonstrated its potential to serve as measurement of how well a S1 feature is able to reflect heterogeneous developments within a landscape. Because S1 features that produce high rates of averages agreement are most likely less sensitive towards field specific developments.

The calculation and field specific visualization of DoT revealed that the majority of depicted fields tended to be ahead of the landscape. Because this study relied on phenological statistics at the state level such bias was expected. Nevertheless, DoT showcases the potential of the landscape-field comparison, as it produces information on growth progress by incorporating multiple S1 features and orbits without the explicit need of in situ data. Because AVA and DoT derive their explanatory power mainly from deviations between field and landscape level time windows instead of in situ observation the impact of biases introduced by point to pixel mismatches inherent to DWD data was mitigated.

When combining DoT with the rate of agreement, flags or labels can be created for fields that exhibit strong negative tendencies and low AVA. Such flags might prove useful for investigating landscapes in terms of their climate vulnerability, resilience or differences in management practice, especially when agents, such as policymakers or consultants who may lack local knowledge are involved. Furthermore, these flags can serve as an initial indicator for pinpointing areas to start ground surveys on crop damage induced by extreme weather events, outbreaks of diseases or pest infestation (Ajadi et al., 2020; Singh et al., 2007; Stobbelaar et al., 2004).

During the explanatory analysis of outliers and deviations across multiple years and crop types, no single, dominant descriptor was identified. The crop specific analysis yielded a slightly stronger emphasis on TWI and OCS/SOC instead of raw elevation. However, as previously mentioned, no dominant descriptor emerged, further supporting the indication of a complex system behind the spatial and temporal occurrence of time windows. Despite this complexity, the model accuracy as indicated by MAE and RMSE (around one to two outlier occurrences) suggests a high prediction potential when combining all environmental descriptors. From a data perspective, it seems logical that these descriptors ranked highly in most explanatory analyses, given that terrain (Yang et al., 2011), soil organic carbon (Santos et al., 2023) and wetness or moisture (Esch, 2018; Pichierri et al., 2018) have well documented impacts on SAR signals. However, their impact on phenological development is highly dependent on the phenological stage and crop type (Canisius et al., 2018).

The GDD-based range of trackability (see chapter 2.6.4) and the initial uptake of signal change provided a different perspective on the accuracy of time series-based monitoring, especially when describing uncertainties. Considering the complexity of S1 time series (with varying viewing geometries, different smoothing algorithms and length of time series) (Harfenmeister et al., 2021b; Löw et al., 2024; Meroni et al., 2021; Qadir et al., 2023), such a measure can be converted to a quality mask at the field level. As mentioned earlier in of the discussion, in situ observations of phenological development are rarely a binary issue. Therefore, a range of uncertainties seems more appropriate for describing the semi-continuous (Meier, 2001) process of crop development, rather than simple temporal differences. However, constant and extensive monitoring of landscapes and their phenological development may improve overall tracking capabilities (d’Andrimont et al., 2022; Tran et al., 2022). Interestingly, the exemplary time line of tracking ranges (see Figure 6) indicate, that drier years produce smaller tracking ranges and thus more accurate predictions. This contrasts with the findings of Schlund and Erasmi (2020), who concluded that tracking quality of InSAR coherence decreases during years of drought. It is however important to note that Schlund and Erasmi (2020) conducted a single orbit study while the present study combined orbits 95 and 146 for winter wheat. Therefore, it is assumed that the effect of different viewing geometries on InSAR coherence result in a larger tracking range within a multi-orbit framework.

Compared to AVA and DoT, the tracking range and trackability are more dependent on the availability of phenological in situ observations. But this study demonstrated that statistics at the level of federal state already provide sufficient information to calibrate the ranges, thanks to the attached GDD baseline and its translation into phenological progress. Hence the point to pixel mismatch within the DWD data is also considered of minor importance.