Satellite-Based Frost Damage Detection in Support of Winter Cover Crops Management: A Case Study on White Mustard

Abstract

Cover crops are grown in order to provide agro-ecological services and must be terminated before planting the subsequent cash crop. Winterkill termination (by frost damage) depends on the interaction between crop frost hardiness, temperatures and the development stage reached at the time of sub-zero temperature exposure. Remotely sensing intensity, timing and spatial variation of cover crop frost damage can be useful for modeling and planning purposes. Therefore, in this study Sentinel-2 vegetation indices were employed in order to detect frost damage in four white mustard (Sinapis alba L.) fields located in Northern Italy. We estimated the starting date of frost events by means of vegetation indices (EVI, NDRE, NDVI, MMSR, and CCCI); we quantified and mapped frost damage at the sub-field level, using ground-based frost damage measurements carried out during the 2021/2022 season. As to frost damage quantification, MMSR outperformed the other VIs followed by CCCI and EVI (R² > 0.55). The adopted procedure to detect starting dates of frost events was successful in most cases, with a one-day and a four-day delay in the two best cases (NDRE). Finally, maps of frost damage were consistent with its observed spatial variation. We demonstrated that it is possible to employ vegetation indices in order to detect cover crop frost damage and thus assessing cover crop winterkill termination efficiency in the field. Further research is needed, involving additional field monitoring of white mustard in more diverse conditions, and extension of the calibration, as well as validation.

1. Introduction

Cover crops are defined as crops grown during bare soil periods between two cash crop growth cycles; these crops are not cultivated to be harvested but rather with the purpose to obtain agro-ecological benefits such as nitrate leaching reduction, weed growth reduction, soil organic matter increase and soil structure stability improvement [1]. Cover crop management normally consists of sowing and termination (i.e., destruction), because in most cases irrigation, fertilisation and weed control are not needed. Termination is carried out before seed bed preparation for the following cash crop. Cover crop management is aimed at maximizing cover crop benefits and minimizing costs, since its biomass is left in the field after termination. Termination can be carried out either mechanically by disking, plowing, or roller-crimping [2], or chemically by herbicide application [3]. Frost-sensitive cover crops, when sown during autumn, may also be efficiently terminated by winter frost damage [4]. Frost termination (winterkill) is particularly convenient in conservation cropping systems (where the choice among mechanical termination methods is limited), and in organic cropping systems (where chemical herbicides are not allowed). The occurrence and the efficiency of winterkill are strongly dependent on crop frost hardiness, environmental temperature, and the development stage reached at the time of sub-zero temperature exposure. Therefore, early or late cover crop planting influences winterkill termination occurrence.

In this context, the detection of cover crop frost damage, by identifying the timing of frost damage events and the affected field area, represents a functional tool to support crop field management. Retrieving the frequency of winterkill termination occurrence, in a specific cover crop cultivation area and considering early vs. late sowing dates, will provide farmers useful information about the necessity to carry out herbicide treatment or employ other mechanical termination methods. Therefore, it will improve planning of the sowing date for the following cash crop in the rotation. Furthermore, it will enable evaluation of cover crop spring growth scenarios (absent in the case of winterkill) and residue management scenarios (eased by early winterkill). Remote sensing techniques represent the most powerful tool to detect and spatialize crop status on large areas. One of their most widespread applications consists in using satellite images in order to characterize crop status [5] through the retrieval and the analysis of vegetation indices [6]. Vegetation indices are algebraic combinations of two or more values of reflectance measured in specific regions of the electromagnetic spectrum. Different vegetation indices are used to detect specific crop stressors [7]; however, none of them can serve as a specific indicator of single stress, because more stressors similarly affect vegetation indices values [8]. The detection of frost damage by vegetation indices has been addressed by different works [9,10,11,12]; they tested and identified different vegetation indices mainly on winter wheat. Some works are available for sugarcane [13,14], coffee [15] and oilseed rape [16]. On the basis of the abovementioned studies, we conclude that several vegetation indices are suitable for frost detection; these are mostly based on red, green and NIR bands. All studies established a relationship between the field-average vegetation index (VI) and the frost damage assessed visually by field surveys or by crop yield measures. Nonetheless, no studies of cover crop frost damage assessment by remote sensing have been published yet, and only a few works have focused on how to provide spatialized information about areas affected by frost damage [13,16,17,18]. Satellite-based vegetation indices were used in the above-cited studies mostly to the purpose of territorial analyses, whereas within-field variability of frost damage was retrieved by unmanned aerial vehicles (UAV) [15], which are more expensive and cover smaller areas compared to satellites. Nonetheless, UAV-derived images have a higher resolution than those derived from the satellite. Finally, only a few attempts have been made to automatically identify time windows of frost events [11,14,19,20]; they were mostly based on the comparison between the selected vegetation index dynamics in a certain year with the dynamics of a reference year without frost damage; these approaches were developed to assess frost effects on the yield of cash crops and none of them specifically focused on cover crops.

Therefore, we employed Sentinel-2 vegetation indices to detect the timing and intensity of frost damage on four fields in Northern Italy cultivated with white mustard (Sinapis alba L.); this species was selected among the most widespread winter cover crops in the area because of its low frost tolerance. We aimed at developing and testing a method to retrieve the onset of frost damage events by vegetation indices, and at quantifying, at field level, the affected areas, by means of spatial analysis of multi-temporal images.

2. Materials and Methods

2.1. Study Area and Measurements

2.1.1. Experimental Sites

The field monitoring campaign was carried out during the 2021–2022 autumn-winter season in commercial farm fields distributed in four sites in Lombardy region (comprised between 45°16′ N and 45°36′ N, and between 9°22′ E and 10°17′ E). Sowing dates ranged from the beginning to the end of September (as reported in

Table 1. Monitoring campaign sites, sowing dates, white mustard cultivars and seeding rates.

Site (Province)	Sowing Date	Cultivar	Seeding Rate (kg/ha)
Alfianello (Brescia)	5 September 2021	Maryna	15
Trezzo sull’Adda (Milano)	15 September 2021	Maryna	18
Fontanella (Bergamo)	15 September 2021	Borowska	15
Ghedi (Brescia)	18 September 2021	Zlata	15

). White mustard seed bed preparation was generally carried out with minimum tillage techniques or with direct seeding; soils in two sites received slurry application a few days before sowing.

2.1.2. Ground Measurement of Frost Damage

Frost damage to aboveground biomass (AGB) tissues was quantified by means of periodic AGB sampling. The monitoring campaign started on 16 November 2021 (before the onset of night frosts) and samplings were repeated every 10–15 days until 11 January 2022. Each sampling date was scheduled to be coupled with one of Sentinel-2 passages over the monitored field. Biomass sampling was performed in five different plots at each site. Frost damage was visually evaluated by detecting frost damage symptoms reported in the literature for other Brassicas [21]: leaf tissues water-soaked, translucent and limp appearance (caused by the disruption of cell walls and membranes due to ice crystals formation, or due to the dehydration of undamaged frozen cells), stem breakage and the consequent tissue senescence. During this evaluation we ruled out other causes that can lead to similar symptoms. No stressors other than frost (i.e., no biotic stress, no water and nitrogen deficiencies) were detected in the monitored fields. Biomass was collected from a 1 m² sampling area for each plot and was then divided in two components: a frost-damaged component (AGB with symptoms) and a non-damaged one. To measure the AGB dry matter content (% DM) of each component, samples were oven-dried until constant weight at 105 °C. In the case of fresh samples weighting more than 1 kg, only a representative sub-sample was dried. A ground-measured frost damage index (GI, Equation (1)) was calculated for each combination of sampling date and monitoring site, by dividing the frost-damaged AGB (t DM ha⁻¹) by the total amount of AGB (t DM ha⁻¹), the latter consisting in the sum of the two sampled AGB components (frost damaged and non-damaged ones).

Ground-measured frost index (GI %) = (frost damaged AGB)/(total AGB)

(1)

Therefore, GI ranges between 1 (all the AGB sampled in the reference area is damaged) and 0 (all the AGB sampled in the reference area is not damaged). The GI was also used to determine the ground-based time window of frost damage (See Section 2.2.1). The environmental temperature was monitored in two sites (Alfianello and Trezzo sull’Adda) by means of three temperature sensors (Watchdog external sensors 3667, Spectrum Technologies, Aurora, IL, USA), placed at 2.0 and 0.5 m above the soil surface to record air temperature, and at 0.1 m depth to record soil temperature. In Ghedi, air temperature at 2 m height was monitored by a weather station next to the field (Vantage Pro2 GroWeather, Davis Instruments, Hayward, CA, USA), while Fontanella data were collected from a weather station belonging to the Regional Agency for Environment Protection (10-km away).

2.2. Satellite-Based Frost Damage Detection

The following software applications were employed: Google Earth Engine (GEE; [22]); QGIS version 3.16 (QGIS Development Team, QGIS.org); R version 4.0.5 ([23], 31 March 2021) and Management Zone Analyst (MZA version 1.0.1, USDA; [24]).

2.2.1. Retrieval of Sentinel-2 Images and Calculation of Vegetation Indices

The Sentinel-2 products (L2A level) used came from the ‘COPERNICUS/S2_SR’ dataset provided by GEE, which was used to elaborate and download the images of the fields. Since cover crops were always sown in September, the time series used in this work refer to the period between 1 September 2021 and 31 January 2022. The actual dates of the images used for frost damage detection are listed in Table S1 of the Supplementary Materials. After applying the time filter, images were processed in three steps: (i) cloud masking, (ii) calculation of vegetation indices, (iii) vectorization of images as data points. Cloud masking detected clouds by applying a filter that selected images with cloud cover under 30%. Then, the QA60 band classified every pixel in each image by cloud presence (0 for no opaque or cirrus clouds, and 1 for the presence of opaque or cirrus clouds); a filter was then applied, and only cloud-free pixels were retained. Five vegetation indices were selected on the basis of previous literature on remote sensing for monitoring frost damage and then tested (

Table 2. Vegetation indices used in this work and their references. The original formulas are reported indicating the reflectance bands (R_band) or the wavelengths (expressed in nm). The formulas using Sentinel-2 bands identify the corresponding satellite band (B) used.

Vegetation Index	Vegetation Index (Name)	Formula	Formula Using Sentinel-2 Bands	References
EVI	Enhanced Vegetation Index	$2.5·\frac{{\mathrm{R}}_{\mathrm{NIR}}-{\mathrm{R}}_{\mathrm{RED}}}{{\mathrm{R}}_{\mathrm{NIR}}+6·{\mathrm{R}}_{\mathrm{RED}}-5·{\mathrm{R}}_{\mathrm{BLUE}}+1}$	$2.5·\frac{\mathrm{B}8-\mathrm{B}4}{\mathrm{B}8+6·\mathrm{B}4-5·\mathrm{B}2+1}$	[11,12]
NDRE	Normalized difference Red-Edge	$\frac{{\mathrm{R}}_{\mathrm{NIR}}-{\mathrm{R}}_{\mathrm{RED}\mathrm{EDGE}}}{{\mathrm{R}}_{\mathrm{NIR}}+{\mathrm{R}}_{\mathrm{RED}\mathrm{EDGE}}}$	$\frac{\mathrm{B}9-\mathrm{B}5}{\mathrm{B}9+\mathrm{B}5}$	[12,15]
NDVI	Normalized Difference vegetation index	$\frac{{\mathrm{R}}_{800}-{\mathrm{R}}_{670}}{{\mathrm{R}}_{800}+{\mathrm{R}}_{670}}$	$\frac{\mathrm{B}8-\mathrm{B}4}{\mathrm{B}8+\mathrm{B}4}$	[10,11,12,15,20]
MMSR	Modification to Modified Simple Ratio	$\frac{800\mathrm{nm}-670\mathrm{nm}-1}{\sqrt{800\mathrm{nm}-670\mathrm{nm}+1}}$	$\frac{\mathrm{B}8-\mathrm{B}4-1}{\sqrt{\mathrm{B}8-\mathrm{B}4+1}}$	[15]
CCCI	Canopy Chlorophyll Content Index	$\frac{\mathrm{NDRE}}{\mathrm{NDVI}}$	-	[12]

).

Finally, a filter for each site was applied by using the shapefile of field boundaries. For each field, the satellite images referring to the time series for the five vegetation indices were translated into a text file (.csv) of points corresponding to the centroids of each pixel (creating a grid of 10 m × 10 m spatial resolution) and downloaded.

2.2.2. Intensity Quantification and Identification of Starting and Ending Dates of Frost Damage Events from Vegetation Indices

Using vegetation indices, we estimated two essential quantities regarding frost damage: firstly, its intensity; secondly, its starting and ending date.

Intensity. For each sampling date, five pairs of GI and VI values were obtained from the five sampling plots in each site. Each VI plot value was retrieved from the pixels belonging to the plot sampled area. Linear regressions were estimated, separately for each VI, between GI and normalized vegetation indices (nVIs). The values of VIs were normalized for each site by dividing them by the average VI of the first sampling date (considered as the baseline VI of a non-damaged crop) [25].

Beginning of frost damage. For the starting date of frost damage, the theoretical basis of the procedure relies on the principle that the curve of the vegetation index can display a drop due to the effects of frost damage on the above ground biomass of the crop [11]; however, other stressors can impact the vegetation index. The stress phase can also be recovered during the season, with a corresponding increase of the vegetation index [11]. To distinguish a decrease of the vegetation index determined by the frost from a decrease due to other causes, the daily minimum temperatures must be recorded and compared to critical thresholds (set for the studied species). When daily minimum temperatures are lower than critical thresholds, the observed decrease of the index is ascribed to a frost damage event. If no increase of the vegetation index is detected after that moment, the start of the frost event is finally set. Based on these considerations, we have developed two methods for identifying the start and the end of a frost event. The first method (‘peak-based method’) comprises the following steps:

(1)

Determining the set of peaks from the curve of the vegetation index according to its first derivative.

(2)

Determining the time windows in which daily minimum temperature (at 2 m height) is under the 0 °C threshold.

(3)

Retaining the dates of the peaks that fall in the time windows determined in step 2.

(4)

Selecting the peak with the highest vegetation index among those identified in step 3. The time window associated with this peak is the period when the frost damage has occurred.

The second method (a semi-supervised procedure called ‘rate of decrease method’) is applied when the ‘peak-based method’ is not able to find frost damage time windows. Since the vegetation index decreases over time, the first derivative calculated at step 1 is taken into account and, if a negative rate (dVI/dt) is observed within a time window of daily minimum temperatures (at 2 m height) under 0 °C, the start of frost event is set on the day in which the first negative dVI/dt was recorded. Then the timeframe of the frost event is closed on the last consecutive day in which a negative daily minimum temperature was recorded.

A linear mixed model was applied to each combination of site and VI (p < 0.05), to assess the statistical relevance of the differences between VIs values measured at the considered dates. Post-hoc tests were performed employing Sidak correction.

In order to measure the accuracy of the proposed procedure in identifying the first frost event causing damage, we ascertained whether the frost time window computed from vegetation indices occurred earlier or later than the frost time window measured on the ground. The ground-measured time window started on the last date when GI was equal to zero and ended on the first date in which GI was higher than zero.

The accuracy of satellite-based detection of the onset of the first frost event was determined by graphical indicators (Figure 1) considering the overlapping between the satellite-based time window and the ground-measured time window. The starting date of the frost event detected from the satellite was considered in advance when the satellite-based time window included the ground-based time window; it was considered late when the satellite-based time window was totally included in the ground-based one (high accuracy) or when the satellite-based time window was partially included and it closed after the first date with a GI higher than 0 (low accuracy).

2.2.3. Spatial and Temporal Variation of Frost Intensity

We aimed to test which combinations of pixels and dates were most similar in terms of frost damage, as represented by vegetation indices. To do so, multi-temporal clustering was applied to each field. A 3D matrix of the field was produced by combining maps referring to different dates, so that the first two dimensions were longitude and latitude, with the third dimension being time. Then, the time series of each pixel was extracted and unfolded into one vector per vegetation index. Therefore, each vector contained the sequence of the vegetation index values for all pixels and all dates. Then, k-means cluster analysis was performed on this vector using MZA. The statistically significant number of clusters was retained when Fuzziness Performance Index (FPI) and Normalized Classification Entropy (NCE) index reached their minimum. In principle, each cluster can contain pixels from different parts of the field and different dates, depending on the similarity of the vegetation index between points in space and moments in time. Then, the dataset was re-folded again to obtain the maps of the resulting clusters. The descriptive statistics were calculated to describe each cluster.

Finally, the linear regressions built between the values of the vegetation indices and the GI were applied to each field in order to map the dynamics of frost damage.

3. Results

3.1. Frost Damage Intensity Quantification from Vegetation Indices

The linear regressions between each normalized vegetation index (nEVI, nMMSR, nNDRE, nCCCI, and nNDVI) and the ground-measured frost index (GI, %) are reported in Figure 2. The simple linear regressions obtained with the indices nMMSR, nCCCI, nEVI and nNDRE showed a coefficient of determination (R², unitless) values ranging between 0.46 and 0.64 and were highly significant (p-values < 0.001). The simple linear regression obtained with NDVI reported the lowest R² (equal to 0.21), but it was highly significant (p-value < 0.001).

These results indicate that all the models consisted in inverse type relationships: GI increased when the VI decreased. The only exception was the site with the earliest sowing date, Alfianello, characterized by a VI decrease that did not correspond to an increase of GI; these details of Alfianello are discussed in Section 4.

3.2. Identification of the Starting and Ending Dates of the Frost Events from Vegetation Indices

First method (‘peak-based’). The starting dates of the frost events identified by means of the ‘peak-based method’ applied with MMSR, CCCI, NDVI, and NDRE are reported in

Table 3. Start and end dates of the frost events for each VI × site combination using the ‘peak-based method’. EVI did not identify frost damage.

Site	Alfianello		Fontanella		Ghedi		Trezzo Sull’adda
Date of minimum air temperature	30 November 2021		30 November 2021		30 November 2021		29 November 2021
Minimum air temperatures at 2.0 m/0.5 m (°C)	−4.3/−5.1		−2.5		−2.8		−1.8/−3.1
	Start	End	Start	End	Start	End	Start	End
MMSR	**		7 January	26 January	29 November *	1 December *	7 January	26 January
NDRE	28 November *	1 December *	3 December	23 December	29 November *	1 December *	3 December	23 December
NDVI	3 December *	23 December	3 December	23 December	10 December	12 December	3 December	23 December
CCCI			3 December	23 December			3 December	23 December

* Start and end dates marked with ‘*’ are comprised in the time window of the first appearance of frost damage (GI > 0). ** For empty cells, the method applied did not identify a frost event.

. The method applied to vegetation index EVI did not allow the identification of the starting date in all sites. The measured air temperatures (at 2 m height) that were used in the procedure are reported in the Supplementary Materials (Figures S1–S4).

Generally, the dates of the frost event retrieved from satellite data were consistent with ground observations: in most cases they indicated time ranges including the date in which the minimum temperature was recorded in each site and after which frost damage symptoms were observed on the ground (with GI values greater than 0%). The minimum temperatures of the period between the sowing date and the detection of the first symptoms (GI > 0) are reported in Table 3.

The comparison of the frost time windows estimated from vegetation indices with the time windows measured on the ground (Figure 3, Figure 4 and Figure 5) indicates, in eight VI × site combinations, a delay of the satellite-based frost damage assessment, ranging from 3 to 39 days, with an average of 17.5 days. On the contrary, in five VI × site combinations satellite-based dates of frost events were found before the date on which the minimum temperature was recorded: the advance ranged from 1 to 2 days, with an average of 1.6 days.

The vegetation index MMSR allowed the identification of frost events in three out of four sites (as reported in Figure 3). In Ghedi the time window of frost damage occurrence was consistent with ground observation: it started 1 day before the date of the minimum temperature of the period between sowing and the appearance of the first symptoms. In Trezzo sull’Adda and Fontanella the time window derived from MMSR was delayed by 39 and 38 days, respectively.

The vegetation index NDRE allowed the identification of frost events in all sites (as reported in Figure 4). For NDRE the average difference between satellite-derived and ground-measured start date of frost event was equal to 0.25 days, and lower than the average differences reported for the other VIs. In Trezzo sull’Adda the difference was the highest (4 days of delay), while in Ghedi it was the lowest (1 day of advance).

Considering the vegetation index CCCI, it identified the frost events in two sites out of four (as reported in Figure 5). The estimated starting of frost events showed a delay more relevant in Fontanella and lower in Trezzo sull’Adda: 4 day and 38 days, respectively.

The results of the application of the ‘peak-based method’ with NDVI are shown in the Supplementary Materials (Figure S5).

The second method (‘rate of decrease’). In several VI × site combinations, the VI was characterized by a slow decrease without peaks (e.g., MMSR in Alfianello, Figure 3). Therefore, the ‘peak-based’ method failed, as illustrated in Figure 3, Figure 4 and Figure 5. For these cases, the second method (‘rate of decrease’) was used to identify the start and end date of frost events. The results are reported in

Table 4. Start and end dates of the first frost events for each VI × site combination obtained using the ‘rate of decrease’ method. Dates and minimum temperatures are reported in Table 3.

Site	Alfianello		Fontanella		Ghedi		Trezzo Sull’adda
	Start	End	Start	End	Start	End	Start	End
EVI	29 November *	1 December *	29 November *	1 December *	29 November *	1 December *	7 December	23 December
MMSR	29 November *	1 December *	29 November *	1 December *	7 December	7 December	7 December	23 December
CCCI	29 November *	1 December *	29 November *	1 December *	29 November *	1 December *	22 December	23 December

* Start and end dates marked with ‘*’ are comprised in the time window of the first appearance of frost damage (GI > 0).

for the VIs that failed to identify the frost events by means of the ‘peak-based’ method in all sites (EVI) or in some of them (MMSR and CCCI).

Trends over time for EVI and its rate of variation are reported in Figure 6, while the results for MMSR and CCCI are in Figures S6 and S7 of the Supplementary Materials. Time windows identified with the ‘rate of decrease’ method were generally shorter than those obtained with the ‘peak-based’ method illustrated before. Furthermore, the frost events obtained with the ‘rate of decrease’ method were more precise compared to the ‘peak-based’ method. The error in the estimation of the start date was, for four VI × site combinations, a delay (ranging from 7 to 23 days, with an average of 11.5 days), while in eight VI × site combinations it was an advance (always equal to 1 day).

EVI is in grey. The starting date of the first frost event estimated with the ‘rate of decrease’ method based on EVI is a vertical light-blue dashed line, while the ending date is a vertical light-blue dotted line.

The results of the application of the ‘rate of decrease’ method with NDVI are shown in the Supplementary Materials (Figure S8).

The significant differences of VIs values between dates, within each VI × site combination, are reported in Table S2 of Supplementary Materials. The VIs values corresponding to the peaks identified by the ‘peak-based method’ were in all cases (except for NDRE peak in Alfianello) significantly higher than the VIs values measured in the following date. Alfianello NDRE peak value, recorded on 29 November, was comparable to the NDRE values measured during the following two dates (7 and 12 December): this VI value decreased significantly only after 22 December. As to EVI, MMSR and NDRE, the VIs values corresponding to the first date with negative dVI/dt identified by the ‘rate of decrease method’ were in all cases (except for MMSR and NDRE, respectively in Ghedi and Alfianello) significantly lower than the VIs values measured during the preceding date. Concerning CCCI trends, this behavior was observed in Trezzo sull’Adda, while for NDVI it was observed in Fontanella and Trezzo sull’Adda.

3.3. Spatial and Temporal Variation of Frost Intensity

The mapped results after multi-temporal clustering of MMSR and CCCI (the VIs with the best correlation with GI) are shown for the field in Ghedi, representative of the behavior of this VI at all sites (Figure 7).

The multi-temporal clustering classified differences of VI mainly among dates and not within the field, indicating that sharp differences of the VI are due to frost damage and not to spatial variability detected by the satellite. The index MMSR clearly showed a discontinuity between 12 and 22 December (Figure 7), when the GI estimated from ground measurements was higher than 60%: on that date, the map of MMSR suggests that the field was mostly damaged (cluster 2). The CCCI was able to identify some spatial differences within the field. Nonetheless, after 22 December it showed a decrease in its values corresponding to 100% of measured GI (Figure 7).

Due to its high sensitivity to frost damage, MMSR was used to study the spatio-temporal distribution of the assessed frost damage in the fields of Ghedi (Figure 8) and Alfianello (Figure 9). The field in Ghedi was chosen as it was representative of the sowing date of common practice, while Alfianello represents an earlier sowing date. The maps of the estimated GI were consistent with the spatial pattern observed on the ground during monitoring. At both sites the frequency of high GI increased over time. As a consequence of sub-zero temperatures recorded at the end of November (Table 3), after 29 November most of the pixels of the two fields showed an estimated GI > 50%. Nonetheless, it is only in January that the two fields recorded 100% of estimated frost damage; however, the field in Alfianello never reached GI of 100% of ground-based frost detection (it reached about 70%, Figure 3). Despite the ground-truth, it should be noted that the field in Alfianello had a higher rate of decreasing GI with most pixels showing GI greater than 66% from December 12 even if low levels GI were measured.

4. Discussion

Vegetation indices and spectral bands. We tested the feasibility of employing various vegetation indices to identify the timing, intensity and spatial variation of cover crop frost damage. The results showed that VIs calculated from reflectance values in the red and red-edge regions were the most correlated to frost damage: CCCI outperformed the other VIs, followed by MMSR and EVI. Unlike other VIs formulas, EVI formula also includes the blue band (Table 2). The fact that VIs based on blue and red bands can quantify the effects of frost damage is in agreement with [10], which reports the ‘purple region’ of the electromagnetic spectrum as suitable to detect frost damage on leaves. The same work indicated the ‘green-edge’, the portion of the spectrum between green and red bands, as the most adequate to detect frost damage on leaves; as Sentinel-2 is not equipped with sensors in that range we could not test this region. Nonetheless, the effects of frost damage on white mustard plants caused a change in color from green to yellowish-green and yellow, substantiating the evidence of MMSR and CCCI success. Other works, at higher spatial resolution, showed similar good results with these red-based VIs compared to green and red-edge bands-based VIs [15]. Nonetheless, the relationships (Figure 2) to predict GI from VIs are still to be improved because they are based on a relatively small dataset. Therefore, more data are needed to validate these regressions to be used in practical applications; however, the significant differences detected (between peaks values and the value on the following date, and between the VI value of the first date with negative dVI/dt and the preceding date) suggest that the two proposed methods are reliable and able to provide significant information. Indeed, these results demonstrate that during the identified frost time windows, statistically relevant decreases of the VIs values were observed.

Vegetation indices and their trends over time; moreover, the collected data showed that good results in the regression analysis were obtained by indices with a smooth trend (MMSR and CCCI), i.e., with no sharp changes in their values over time (Figure 3 and Figure 5). Due to the good results in the regression analysis, VIs were used firstly to identify the date of the start of frost events, and secondly to map frost damage.

The ‘peak-based method’ to identify the start date applied to MMSR, CCCI, and EVI did not lead to good results, while it was successful with NDVI and NDRE, which were characterized by more spikes (Table 3). Indeed, according to the literature, NDVI and NDRE are more sensitive to frost damage compared to other indices, with a sharper decrease in their values following the damage [14,16,18,26]. When EVI was tested in previous studies, the spatial resolution of satellite images was lower than it is nowadays [16]. In this work, the adopted spatial resolution was 10 m and the rate of decrease method with EVI and MMSR allowed the identification of frost events starting dates with delays under 23 days, even if these VIs had more smoothed dynamics (Figure 3 and Figure 6) compared to the other VIs; these results open the opportunity of using other VIs than NDVI to successfully detect crop frost damage; however, it must be emphasized that all the VIs in Alfianello decreased after mid-December, even if this event did not correspond to an increase in the ground-measured GI (Figure 3); it was partly due to the fact that the plants in the south part of this field, on 29 November, were at the flowering stage with fully opened flowers. Here, the yellow petals of the flowers caused a decrease in the values of the VIs, already seen in the literature [27], which was erroneously ascribed to frost damage in our regressions; this issue cannot be overcome by the use of VIs alone, requiring phenological models to estimate the crop development stage.

Spatial variation of frost damage. We tested VIs ability to identify frost damage spatially at the sub-field level, as there was a lack of dedicated studies and only few works are reported in the literature [15]. Thanks to the significant regression models calibrated with satisfactory coefficients of determination for some VIs, the results of the spatialization of frost damage (Figure 8 and Figure 9) were encouraging because they reflected the spatial variability observed during the field campaigns. Nonetheless, it must be remarked that the multi-temporal cluster analysis adopted to spatialize the frost damage at the sub-field level showed more differences among dates than within the fields (Figure 7), suggesting that the main driver of the values of VIs was frost damage, not local soil or atmospheric conditions.

Perspectives. Overall, the results obtained showed that the use of satellite imagery for sub-field detection of frost damage on white mustard is possible and can be improved by accounting for crop development stage. Once operational, this method could have an agronomic application in supporting territorial-scale stakeholders (e.g., extension services, farmer associations) to identify winter-sensitive cover crop species on large and distant areas. In the case of white mustard in Northern Italy, early sowing dates (September), which normally lead to winterkill termination as in this study, simplify crop residue management because cover crop biomass in spring is already partly decomposed, and allow early sowing of the following cash crop (end of March-beginning of April) thanks to simplified seed bed preparation. On the contrary, late sowing dates (October or even November) can be chosen when the sowing date of the following cash crop is delayed (e.g., soybean) and, therefore, cover crop growth can continue until mid/end of April; this also avoids that cover crop early winterkill leaves the soil bare during a period with a high risk of leaching and weed growth.

The presented methods could be associated to other tools able to monitor catch crops [28] to further support the stakeholders involved in the Common Agricultural Policy Greening measures definition and application.

Future perspectives involve the validation of the proposed method on a larger dataset, studying its applicability on other cover crop species, and the use of crop models to account for the effect of the development stage on VIs.

5. Conclusions

The application domain of this monitoring study consists in providing a tool to monitor cover crops’ frost damage. We demonstrated that it is possible to employ vegetation indices to detect the onset of white mustard cover crop frost damage, quantify it and map its spatial distribution. Firstly, the ability of different vegetation indices to quantify frost damage in terms of dead biomass was tested. The best vegetation indices were the ones based on red, red-edge, and near-infrared regions of the reflectance spectrum of crop canopy (MMSR, EVI, and CCCI). The results were consistent with the literature and allowed the further step of using vegetation indices to identify the starting dates of frost events. Two methods were used: the analysis of peaks of vegetation indices (‘peak-based method’), and the study of the rate of decrease over time of vegetation indices (‘rate of decrease method’). The two methods, developed in previous works, were simplified by us with good results. In fact, previous literature focused on cereal crops and more complex methods of detection, resulting in difficult application in the common practice. Our results were encouraging and supported by the ground-based detection of frost damage. Finally, the study of the within-field spatial variability of frost damage was performed by using the regression models built in the first place and by attempting a downscaling of previous works, which were applied on a territorial scale. Despite the promising results, the effects of the crop development stage were evident in the only field that reached the full flowering stage, where the areas of the field with white mustard yellow flowers were classified as damaged by frost. Further studies should be carried out by coupling vegetation indices and crop simulation modeling in order to account for the crop development stage, thus leading to a more robust mapping of crop canopies frost damage; moreover, further research should involve additional field monitoring campaigns to improve the quality of the regressions between vegetation indices and ground-measured frost damage, taking into account a wider range of crop development stages, and to validate the proposed frost detection methods.