Near Real-Time Freeze Detection over Agricultural Plots Using Sentinel-1 Data

Abstract

Short-term freeze/thaw cycles, which mostly occur in the northern hemisphere across the majority of land surfaces, are reported to cause severe economic losses over broad areas of Europe and North America. Therefore, in order to assess the extent of frost damage in the agricultural sector, the objective of this study is to build an operational approach capable of detecting frozen plots at the plot scale in a near real-time scenario using Sentinel-1 (S1) data. C-band synthetic aperture radar (SAR) data show high potential for the detection of freeze/thaw surface states due to the significant alterations to the dielectric properties of the soil, which are distinctly observable in the backscattered signal. In this study, we propose an approach that relies on change detection in the high-resolution Sentinel-1 C-band SAR backscattered coefficients, to determine surface states at the plot scale as either frozen or unfrozen. A threshold analysis is first performed in order to determine the best thresholds for three distinct land cover classes, and for each polarization mode (VH, and VV). S-1 SAR data are then used to detect a plot’s surface state as either unfrozen, mild-to-moderately frozen or severely frozen. A temperature-based filter has also been applied at the end of the detection chain to eliminate false detections in the freezing detection algorithm due mainly to rainfall, irrigation, tillage, or signal noise. Our approach has been tested over two study sites in France, and the output results, using either VH or VV, compared qualitatively well with both in situ air temperature data and soil temperature data provided by ERA5-Land. Overall, our algorithm was able to detect all freezing episodes over the analyzed S-1 SAR time series, and with no false detections. Moreover, given the high-resolution aspect of S-1 SAR data, our algorithm is capable of mapping the local variation of freezing episodes at plot scale. This is in contrast with previous products that only offer coarser results across larger areas (low spatial resolution).

1. Introduction

Frozen grounds are soil where part or all of the pore water is found as ice. Frozen surfaces are mostly found in the northern hemisphere, and make up the majority of land surfaces [1,2]. In permafrost, the ground remains completely frozen for at least two consecutive years [3]. Seasonal frozen regions have frost persisting for several months. In contrast, short term freeze/thaw (FT) cycles only last from hours to a couple of days. Permafrost make up to 24% of total land surfaces, and 55% of exposed land surfaces undergo short term FT cycles each year [4]. The FT cycles are sensitive to the fluctuation of the isothermal temperature threshold (0°C) at ground level, and have significant impacts on the climate system, surface runoff, water balance, carbon cycle, and crop growth [2,4,5]. The relationship between FT cycles and crop yield are well understood and have been extensively studied [6,7,8,9,10,11]. FT cycles have an impact on both the physical and biological aspects of the soil in which the plants grow. For example, freezing and thawing causes expansion and contraction that push plant roots up and sometimes out of soil, which damages them. Moreover, frost damage occurs to all kinds of crops, like cereal grains, flowers and seeds, canola pods, legume pods and seeds, flowers, and vines, and can damage whole plant heads in the case of soil freeze, which causes the stems to also freeze [12]. To lower the effects of frost damage, generally, producers often plant crops in a manner so that plant flowers, leaves and grains develop during warmer temperatures. However, due to climate change, warmer temperatures are causing premature vegetation growth, and as a result, crops will more likely be affected by late frost which occurs during spring. In the spring of 2020, a swathe of freezing temperatures overtook much of Europe, causing severe damage to already flowering crops throughout Northern Italy, much of Poland, and the Ukraine. In 2018, crop damage from spring freezing events has been reported as causing severe economic losses over broad areas in Europe and amounting to hundreds of millions of euros [13]. In 2007, the spring frost event in south-central and south-eastern USA caused agricultural damage in the order of 1.6 billion euros [14]. In order to assess the extent of frost damage, and subsequently the potential loss of yield in agricultural areas, FT maps are therefore required. However, such mapping requires long time series of data, with periodic measurements over large geographical areas, which is impractical to collect in situ.

To monitor large areas for FT events, in real-time or near real-time, remotely sensed satellite data are primordial. Satellite microwave radars, radiometers and scatterometers are well suited for such monitoring, and they have been extensively used in the past for the monitoring of surface state parameters, such as freeze/thaw [15,16,17,18,19,20]. Radars are active sensors that illuminate a target area and measure the backscattered energy, while radiometers are passive sensors that record natural microwave emissions. The potential of radars, radiometers and scatterometers to detect surface FT is very promising due to the dielectric properties of the soil, which are very sensitive to the amount and state of water contained in the soil. Therefore, any variation in the state or amount of water (e.g., in the case of freezing) causes significant alterations to the dielectric properties of the soil. These alterations are distinctly observable in the radar’s backscattered signals [21,22,23,24]. With active microwave sensors, a decrease in the order of 3 dB in the backscattered signal can be observed between frozen and unfrozen surfaces over the same area [20,22,23]. For passive microwave sensors, a 10 °C decrease in clay or silt soil temperatures was observed to increase the sensed temperature brightness by 14 (clay) and 28 °K (silt) [25].

In recent years, we have observed strong technical progress and an increase in the number of proposals for sensors in orbit. This is particularly true for synthetic aperture radars (SAR) like the recently launched Sentinel-1 constellation of European Copernicus program, and the arrival of long time series at low resolution from radiometer and scatterometer sensors like SMOS, SMAP, and ASCAT. ASCAT and SMOS have already been proven successful in a large body of studies for the detection of freeze/thaw events [15,16,17,18,26], and are capable of monitoring such events due their high repeat cycles (1–3 days for both sensors), albeit at low spatial resolutions (25 km, and 35 to 50 km for ASCAT and SMOS, respectively). Moreover, SMAP is the first satellite to provide a global freeze/thaw product at low resolutions (9 km) with high repeat cycles (two to three days). The product is based on a temporal change detection approach that assumes that the high variation in the temporal dynamics of the normalized polarization ratio (NPR) of the brightness temperature is related to a large change in the dielectric constant (passing from frozen to unfrozen condition, and vice-versa). However, due to the low resolution of these products, they are unable to detect freeze/thaw events at plot scales, which is the required scale for agricultural surfaces. Currently, C-band data provided by Sentinel-1’s (S1) two-satellite constellation offer a good opportunity for the operational mapping of frozen agricultural areas due to S1’s free and open access data at high spatial resolutions (10 m × 10 m) and high revisit times (6 days or less over Europe).

The algorithms developed for the detection of frozen and unfrozen surface states are based either on passive or active microwave sensors. Passive sensor-based algorithms can be classified into four main classes. These algorithms include: (1) dual-indices [25,27], which use thresholds over temperature brightness and spectral gradient parameters for the detection of freeze/thaw events. (2) New dual-indices [19] use thresholds on the standard deviation of horizontally polarized brightness temperatures, and vertically polarized brightness temperature. (3) Decision trees [28,29], using brightness temperatures. Finally, (4) the Discriminant Function Algorithm (DFA) is a recently developed method by Zhao et al. [30], which relies on temperatures brightness and quasi-emissivity. DFA has also recently been shown to have a better freeze/thaw detection accuracy than the other three methods according to Shao et al. [31]. Freeze/thaw detection algorithms using active microwave sensors rely on change detection and can also be classified into three main classes: (1) threshold algorithms [16,20,32,33,34], where the frozen/unfrozen states are flagged using a single backscattered value compared to a predefined threshold; (2) moving average algorithms [35,36,37,38], which flag a frozen/unfrozen surface state based on the difference between the backscattered value of the current day and the mean backscattered values of n previous days; and finally, (3) frozen/unfrozen surface detection using probabilistic models [15,17], based on a variant of hidden Markov models to compute the probability of snow, and FT events using K_u – and C-band scatterometer data.

Recently, Baghdadi et al. [20] used a thresholding algorithm on the freely available C-band synthetic aperture radar (SAR) Sentinel-1 data, with a constant threshold of 3 dB, to successfully map frozen/unfrozen surfaces. The algorithm calculated the difference of the backscattered S1 signal (σ°) acquired on a specific date and its immediate neighbor not considered as frozen. This approach presented some limitations, where some false negatives (frozen surface detected as unfrozen) were detected when, for instance, the SAR returned signal over wet snow (temperature above 0 °C) will already be 3 dB lower in comparison to the returned signal over unfrozen surfaces. Thus, in the event of a post-freeze episode, the freeze event will not be captured. The algorithm also considered a stable soil roughness over time. A variable roughness could also produce false negatives, as high roughness will increase the return signal [39], and the difference Δσ° (σ° with freezing and high roughness −σ° with thawing and low roughness) could be well lower than the 3 dB threshold. Finally, the study was carried out with data during the winter season, and over surfaces with lower vegetation cover; therefore, it would be interesting to develop an algorithm that works over different vegetation types.

In this study, we propose a refined algorithm based on Baghdadi et al. [20] that is capable of detecting frozen/unfrozen episodes in near real-time, with a temporal resolution of at least once each six days, at plot scale, and over different vegetation types. To overcome some of the limitations encountered during the previous study, two additional datasets will be used: (1) in situ temperature data using the measured temperatures from a weather station network for the refining of the detection results and (2) land cover classes. The inclusion of temperature data will allow the removal of some false positives (unfrozen surface detected as frozen) that could occur in spring, due to specific events such as irrigation [40,41], plowing [39] or vegetation growth [42]. In effect, as the mapping of a surface state as either frozen/unfrozen on a given date ‘$d_i$’ is most often based on the surface state from previous dates, a difference between the radar signal at date ‘$d_i$’ and the radar signal acquired previously can be equal or less than −3 dB. This difference is observable if at previous radar acquisitions, the soil was wet (due to rainfall or an irrigation event), while the radar backscattering coefficient at date ‘$d_i$’ was acquired when the soil was less humid (radar signal increases with humidity [43]). Moreover, the difference between a radar signal at date ‘$d_i$’ and those acquired previously could be less than −3dB if the radar backscatter from previous dates was acquired over rough soil, while at date ‘$d_i$’ the soil was smoother (radar signal increases with roughness [39]). Finally, the difference between a radar signal acquired at date ‘$d_i$’ and backscatter values acquired previously could also be less than -3 dB due to the vegetation growth cycle of some crop types, such as, for example, in the case of wheat [42]. Nasrallah et al. [42] showed that the C-band SAR backscattering signal in VV polarization decreases gradually between the germination phase of wheat, which usually occurs by the beginning of January, and the heading phase occurring between mid-March and mid-April. In the heading phase, the C-band SAR backscattering signal attains extremely low values due to extreme vegetation attenuation. Meanwhile, land cover classes will be used since the SAR backscatter does not decrease during freezing events with the same magnitude over different vegetation types [44,45].

Finally, the generated plots’ surface state maps, which describe a surface state as either frozen or unfrozen, will be compared to the fifth Generation ReAnalysis (ERA5-Land) hourly soil and air temperature data products provided by the European Center Medium range Weather Forecasts (ECMWF) [46]. Several previous studies have used soil and air temperatures as proxies for frozen/unfrozen states [16,17,28,47]. In addition, ERA5-Land predecessor ERA interim has shown close agreement with in situ temperature data [48,49,50]. In addition, the comparison to ERA5-Land’s hourly soil and air temperatures will be performed in a qualitative manner, due to (1) the uncertainties of ERA5-Land’s soil and air temperatures; (2) the exceedingly rare availability of direct ground measurements; (3) the high difference between ERA5-Land’s resolution (9 km), and our study sites’ plot surface areas (several hundredm²), which makes a direct comparison difficult. However, as the aim of this study is to evaluate the potential of S1 SAR data product to detect in near real-time frozen/unfrozen events, comparison to ERA5-Land can identify areas where our algorithm does not perform well. The study sites and the used datasets are described in Section 2, followed by a detailed description of our algorithm for mapping frozen areas at the plot scale in Section 3. The main findings and the discussion are presented in Section 4 and Section 5, respectively. Finally, Section 6 presents the main conclusions.

2. Study Sites and Datasets

2.1. Study Sites

We analyzed around 9 months of data, from 1 September 2018 to 31 May 2019 over two selected sites located in metropolitan France (Figure 1). The criteria of choice for these two sites was that they are agricultural areas, and are likely to be affected by freeze episodes. The first study site (48.86°N, 3.17°E; Site A, Figure 1 left) is situated near Paris and is comprised of 1001 agricultural parcels, with a surface area of 8943 ha, and mean parcel surface of 8.9 ha. In this study site, 80% of the parcels contain cereal grains, 17% meadows, and the rest are orchards and vineyards (roughly 3%). The second study area (48.33°N, 7.43°E; Site B, Figure 1 right) is situated near Strasbourg and is comprised of 11,578 agricultural parcels, with a mean parcel surface area of 1 ha. Here, 68% of parcels contain orchards and vineyards, 23.5% cereal grains, and the rest are meadows. Figure 2 shows the temperature profiles over the two study sites for the period between 1 September 2018 and 31 May 2019. The temperature profiles were acquired from the nearest weather station network to each site, which are ~3 km and ~17 km away from study sites A and B, respectively. Study site A (Figure 2a) shows seven freezing episodes when the temperatures fall below 0 °C and reach a minimum of −3 °C. Six of these freezing episodes occur from late fall (29 October 2018) to the end of the winter season (11 February 2019), and one freezing episode in the middle of the spring season (12 April 2019). Study site B (Figure 2b) shows eleven freezing episodes (temperatures < 0 °C and reach a minimum of −6 °C) in the period from 13 November 2018 to 15 February 15 2019.

2.2. Datasets

2.2.1. SAR Sentinel-1 Time Series

Sentinel-1 is a constellation of two polar-orbiting SAR satellites (S1A and S1B) operating at C-band (ƒ = 5.406 GHz), which were launched in April 2014 (S1A) and April 2016 (S1B). Sentinel-1 products are available free of charge, and can be downloaded one day after their acquisition from the Copernicus website (https://scihub.copernicus.eu/dhus). Over our two study sites, two Sentinel-1 time series data were used. The first S1 time series, from 1 September 2017 to 31 May 2018, was used for threshold calculation. This S1 time series contains 288 S1 images acquired from the Sentinel-1A (S1A) and Sentinel-1B (S1B) satellite constellation in both ascending (afternoon ~18h00 UTC) and descending modes (morning ~06h00 UTC). Of these 288 images, 133 images were acquired over study site A (88 images in ascending mode, and 45 descending mode), and 154 images acquired over study site B (87 images in ascending mode, and 67 in descending mode). The second S1 time series, from 1 September 2018 until 31 May 2019, was used for the mapping of freezing events. The second S1 time series (1 September 2018 and 31 May 2019) contains 263 S1 images acquired from S1A and S1B in both ascending (~18h00 UTC) and descending modes (~06h00 UTC). In the second time series, 130 images were acquired over study site A, of which 43 images acquired in descending mode and 87 images in ascending mode. For study site B, 133 images were acquired, of which 66 images and 67 images were acquired in descending and ascending modes, respectively. Figure 3 shows the repeat cycle of the S1 SAR in ascending “A” and descending “D” acquisition modes over both study sites for the month of October 2018.

For study site A, two images were acquired in the ascending mode at incidence angles of 41.9° and 32.8°, respectively, with a time difference of 2 days, followed after 12 h by an image in the descending mode at an incidence angle of 37.8° (Figure 3a). The pattern was repeated every 6 days for a total of 16 images. For study site B, two images in the ascending mode were acquired in the span of 24 h at incidence angles of 42.1° and 32.8°, respectively, followed by one image in the ascending mode after 6 days at an incidence angle of 32.8°. This pattern was repeated every 12 days (Figure 3b). The acquisition pattern for the images in the descending mode was two images (24 h difference) at incidence angles of 36.6° and 45.2°, respectively, followed after 5 days by an image with an incidence angle of 36.6° (Figure 3b). This pattern was also repeated every 12 days for a total of 16 images in both acquisition modes. The time difference between the acquisition of images in ascending and descending modes varied between two and three days. The entirety of the S1 SAR images database was acquired in the Interferometric Wide swath mode (IW) in both VV and VH polarizations. These images were generated from the high-resolution Level-1 ground range detected (GRD) product, which were acquired at a spatial resolution of 20 m × 22 m (range × azimuth), and then resampled to produce the final product at 10 m × 10 m pixel spacing. The images were calibrated using Sentinel SNAP toolbox, developed by the European Space Agency (ESA). The calibration aimed to convert the image pixel values into backscattering coefficients (σ°) in linear units, while the geometric correction was used to ortho-rectify the SAR images using the 30 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM). The radiometric accuracy of the S1 SAR backscattering coefficients is approximately 0.70 dB (3σ) for the VV polarization and 1.0 dB (3σ) for the VH polarization [51,52]. Finally, the backscattering coefficients ${\sigma }_{\theta }^0$ acquired at different incidence angles θ were normalized to a common reference incidence angle set at θ_ref = 40° using the square cosine correction equation [53,54]:

$${\sigma }_{ref}^0=\frac{{\sigma }_{\theta }^{0}co{s}^2\left({\theta }_{ref}\right)}{co{s}^2\left(\theta \right)}$$

(1)

2.2.2. Land Cover Map

The behavior of SAR backscatter profiles after a freeze event can vary depending on the soil properties and vegetation types [44,45]. Therefore, it can be beneficial for the mapping of frozen/unfrozen areas to study the C-band backscatters dependently of the main vegetation types. The graphical parcel registry (Registre Parcellaire Graphique; RPG) is a graphical declaration system used by farmers in order to provide a fine and annual geo-localized representation of the occupation of agricultural lands in France. The RPG includes around 6 million small parcels (corresponding to approximately 26 million hectares). For this study, we used the 2018 version of the RPG, henceforth named RPG-18, which contains 28 classes of land cover. The classes corresponding to agricultural activities (summer and winter crops, vines, fruit trees and meadows) were used in this study and grouped into three higher-level classes named: cereal grains, orchards and vineyards, and meadows.

2.2.3. Temperature Data

In situ temperatures acquired from the nearest weather station network to each study site (approximately 3 and 17 km from the study sites A and B, respectively) were used to help derive the thresholds for the freezing detection algorithm. For each S1 SAR image, over the two study sites, we matched air temperatures from the time of the S1 SAR acquisition and up to three hours before the acquisition. Then, these temperatures were averaged because freeze/thaw events could have happened at the same time, or a few hours before the image acquisition.

In addition, the ERA5-Land dataset, provided by the European Centre for Medium Range Weather Forecasts (ECMWF) was also used. ERA5-Land is obtained by the modelling of in situ and satellite observations, and provides a large number of meteorological and oceanographic variables [36]. ERA5-Land data are available at a resolution of 9 km globally, in hourly intervals, and it is produced in near real time (three months after acquisition), with quality-assured monthly updates, and preliminary daily updates published within five days of acquisition. For the validation of our approach of frozen area mapping, we used the soil temperature and the 2 m air temperature. The ERA5-Land soil temperature products are composed of four layers, and go to 289 cm in depth. In this study, we used the upmost layer, which goes up to 7 cm in depth (C-band microwaves can penetrate the soil to only a few centimeters). The spatial matching between our plots and the ERA5-Land grid has been achieved by a nearest neighbor approach.

3. Materials and Methods

C-band Sentinel-1 SAR, as all active microwave systems, is sensitive to different physical phenomena that occur in a given landscape. Over bare soils and sparsely vegetated fields, the parameter that governs the C-band backscatter value, in addition to soil roughness, is the soil water content due to its high dielectric constant at microwave frequencies [43,55]. An increase in the soil’s water content will increase the SAR backscatter value [40,41]. However, when ice begins to form in the soil, the dielectric constants decreases to an extent which causes the C-band backscatter values to drop significantly [20,23,24,32]. The backscatter response can also vary significantly with respect to soil and crop types [44,45]. SAR backscatter can increase or decrease between two satellite acquisition dates due to the development of the vegetation cover [42,56]. SAR backscatter is also sensitive to soil roughness, which can increase significantly due to field plowing. It is has been proven that an increase in soil roughness can also increase the backscatter [39,57]. Based on these sensitivities, we propose a surface state detection algorithm that discriminates between frozen and unfrozen surface states. The algorithm is a variant of the threshold algorithm [16,20,32,33,34], combined with a variant of the moving average algorithm, which flags a surface as frozen or unfrozen based on the average of the previously acquired backscatter values [35,36,37,38]. The algorithm developed in the present study is based on the difference between the S1 backscattering coefficient at current date ‘$d_i$’ and a reference S1 backscattering coefficient calculated from S1 backscattering coefficients acquired in the n-days preceding the S1 acquisition at current date ‘$d_i$’. If the difference between the two backscatters (${\sigma }^0$ at date $d_i$ - the reference ${\sigma }^0$) is greater than a threshold that remains to be determined, then it is highly probable that the surface state is frozen at acquisition date ‘$d_i$’ of the S1 SAR backscatter. Moreover, since in the event of freezing, the magnitude of decrease in S1 SAR backscatter is highly correlated with the severity of the freezing event [16], in this study, we used two thresholds to distinguish between mild-to-moderate and severe freezing events. The threshold values were based on the land cover type, as in the event of freezing the decreases in the S1 SAR backscatter are also dependent on the land cover type [44,45]. Our approach can be divided into two independent parts. The first part allows the calculation of the threshold values for each land cover type, and needs to be performed only once, as the same thresholds can be used for different study sites with the same land cover types. The second part consists of the algorithm for the near real-time frozen/unfrozen detection at plot scale, the output of which will indicate the plot state as one of three possible outcomes: (1) unfrozen, (2) mild-to-moderately frozen, and (3) severely frozen.

3.1. Plot Surfaces Frozen/Unforzen Detection Algorithm

To flag a plot (p) as frozen at any given time, the difference between the S1 backscatter over that plot at current date ‘$d_i$’ and a reference backscatter obtained from previous acquisitions needs to be compared against a threshold. In some studies, the acquired backscatter value at a certain date is compared directly to a threshold [16]. Over vegetated areas, C-band SAR backscatter σ° is the contribution of both the soil and vegetation components [55]. The contribution of each one of these two components to the total SAR backscatter differs from one crop type to another. Thus, in the case of freezing, the SAR backscatter of different plots, even for the same land cover, will not decrease with the same magnitude as the characteristics of the soil and vegetation are not exactly the same.

First, a maximum ${\sigma }^0$ is calculated over each plot p from S1 acquisitions acquired n days before, and up to the current date $d_i$ ($max\_{\sigma }^0\left(d_n,p\right)$). The $max\_{\sigma }^0\left(d_n,p\right)$ should not correspond to a freezing event previously detected. The number of days (n) on which the maximum calculation $max\_{\sigma }^0\left(d_n,p\right)$ is carried out should not be too long to minimize the effect of vegetation growth. Therefore, for each plot p at date $d_i$, $max\_{\sigma }^0\left(d_n,p\right)$ is calculated based on the maximum of the ${\sigma }^0$-values corresponding to S1 images acquired not more than 15 days before $d_i$ (at least three images). At date $d_{i+1}$, a new $max\_{\sigma }^0\left(d_n,p\right)$ is calculated in plot p if, and only if, the time difference between $d_n$ and $d_{i+1}$ is greater than 15 days. If not, we use the last calculated $max\_{\sigma }^0\left(d_n,p\right)$. Since $max\_{\sigma }^0\left(d_n,p\right)$ is calculated using the maximum value of backscatters from n previous dates, sometimes, especially after heavy rainfall, or an irrigation event, the backscatter value at date $d_n$ increases abruptly. Hence, to reduce the effects of such influences, when detecting the surface state of a plot p at $d_i$, instead of calculating the difference between the backscatter value at date $d_i$ and only the last calculated $max\_{\sigma }^0\left(d_n,p\right)$, we will instead calculate the difference between ${\sigma }^0$ at date $d_i$, and the mean of the last three calculated $max\_{\sigma }^0\left(d_n,p\right)$. Therefore, after selecting $max\_{\sigma }^0\left(d_n,p\right)$ for plot p at date ‘$d_i$’ _i, we first calculate the mean of the last three calculated $max\_{\sigma }^0\left(d_n,p\right)$. This mean is the reference σ° ($ref\_{\sigma }^0\left(p\right)$) that will be used when calculating the difference between $ref\_{\sigma }^0\left(p\right)$ and ${\sigma }^0\left(d_i,p\right)$ ($\mathsf{\Delta }{\sigma }^0\left(d_i,p\right)$). Finally, $\mathsf{\Delta }{\sigma }^0\left(d_i,p\right)$ will be compared to a threshold range which depends on the land cover type (LC) of plot p and the S1 polarization $\left(\left[A_{lc},B_{lc}\right]\right)$. If $\mathsf{\Delta }{\sigma }^0\left(d_i,p\right)<A_{lc}$, plot p is considered as unfrozen at date ‘$d_i$’. If $\mathsf{\Delta }{\sigma }^0\left(d_i,p\right)\in \left[A_{lc},B_{lc}\right[$, plot p is considered as having a mild-to-moderate freeze occurrence. Finally, if $\mathsf{\Delta }{\sigma }^0\left(d_i,p\right)\ge B_{lc}$, plot p is considered as having a high freeze occurrence.

3.2. False Detection Filter

Due to its dielectric properties, water has a very important impact on the soil’s SAR backscattered signal, as an increase in water content will result in an increase in the backscatter value [40,41,58]. Therefore, irrigation events and rainfall usually cause sudden spikes in the backscatter values. Moreover, the spikes in the backscatter signal can also appear as a result of the increase in surface roughness due to the plowing of agricultural parcels [39]. A false freezing detection could thus be observed following specific events (e.g., irrigation/rainfall or plowing). Assuming a plot state at date ‘$d_i$’, that was preceded at date ‘$d_j$’ by such an event (e.g., irrigation/rainfall or plowing), where ${\sigma }_0$ at date $d_j$ is chosen as the maximum ${\sigma }_0$ $\left(max\_{\sigma }^0\left(d_j\right)\right)$. In such a case, and even though we consider the mean of the last three $max\_{\sigma }^0\left(d_j\right)$ $\left(ref\_{\sigma }^0\right)$, the difference between $ref\_{\sigma }^0$ and ${\sigma }_0$ acquired at date $d_i$ could be of a magnitude that causes the algorithm to detect the plot state at date $d_i$ as being frozen, even if it is not. $max\_{\sigma }^0\left(d_j\right)$ could be high in comparison to ${\sigma }^0\left(d_i\right)$ due to irrigation or rainfall events but also due to plowing. In addition, ${\sigma }^0\left(d_i\right)$ could also be much lower in comparison to ${\sigma }^0\left(d_j\right)$ due to a decrease in soil water content or a decrease in soil roughness (e.g., preparation of soil for sowing). In addition, the difference between $ref\_{\sigma }^0$ and ${\sigma }^0\left(d_i,p\right)$ could be greater than +3 dB if at date $d_i$ there was an important decrease in the radar signal due to the growth cycle of certain vegetation types such as wheat [42].

To minimize such false detections, a filter is applied at the end of the freeze detection algorithm using a temperature threshold from in situ measurements. If a given plot at a certain date $d_k$ is detected as frozen, and the air temperature at the same date is higher than a certain value, the plot state will be set to unfrozen. A temperature of 0 °C was used as a temperature threshold to distinguish between unfrozen and frozen conditions in some studies [59]. Other studies used a more conservative threshold to improve the likelihood of the temperature to represent unfrozen/frozen conditions (e.g., 3 °C threshold in Kim et al. [60]). In the present study, similarly to Kim et al. [60], a threshold of 3 °C will be used to filter out the false detections.

3.3. Threshold Analysis over ${\sigma }^0$

In order to detect surface states over a given plot p at date ‘$d_i$’, we need to calculate a threshold over ${\mathsf{\Delta }\mathsf{\sigma }}^0\left({\mathrm{d}}_{\mathrm{i}},\mathrm{p}\right)\equiv \left(\mathrm{ref}\_{\mathsf{\sigma }}^0-{\mathsf{\sigma }}^0\left({\mathrm{d}}_{\mathrm{i}},\mathrm{p}\right)\right)$ which allows the plot surface to be flagged as being unfrozen, mild-to-moderately frozen, or severely frozen. Since the behavior of the backscatter varies according to the characteristics of the soil and vegetation, a threshold was determined for each land cover class. In this study, crop types are classified into three main land cover classes: cereal grains (LC1), meadows (LC2), and orchards and vineyards (LC3). The thresholds were determined using an approach based on the use of an additional S1 SAR time-series from 1 September 2017, till 31 May 2018. First, the approach consists of calculating the difference ${\mathsf{\Delta }\mathsf{\sigma }}^0\left({\mathrm{d}}_{\mathrm{i}},\mathrm{p},\mathrm{t}\right)$ for each plot p of our three land cover classes, at date $d_i$ when the temperatures t are below 0 °C (temperatures with freezing potential). The temperatures t corresponding to the S1 acquisition dates are acquired from the nearest weather stations. Thus, for each of our three land cover classes (LC) we have a set (S_LC) of ${\mathsf{\Delta }\mathsf{\sigma }}^0$ containing N plots and M dates (when t <0 °C). Then, from each set S_LC we calculate the distribution of ${\mathsf{\Delta }\mathsf{\sigma }}^0$ at each S1 polarization (VH and VV) separately for two temperature ranges, t ∈ [−3 °C, 0 °C] and t < −3 °C. The obtained distributions of each LC at each polarization will be fitted with a normal distribution of mean μ(LC_k,pol_p) and standard deviation σ(LC_k,pol_p), where k is the land cover class (between 1, and 3), and p the polarization, with pol₁ designating VH polarization and pol₂ the VV polarization. In Baghdadi et al. [20], it was shown by simulation that the difference between the backscattering coefficient of unfrozen soil and the backscattering coefficient of frozen soil with soil temperature (t) lower than 0 °C abruptly increases for t between 0 and −3 °C. The radar signal decreased from roughly 3.3 dB for silty clay soil to ~5dB for sandy loam soil when the temperature decreased from 0 to −3 °C. For t between −3 and −20 °C, the radar signal slows its decline to less than 1 dB.. Therefore, in this study, a temperature of −3 °C is thus considered as the temperature at which it is highly probable to have severe freezing. The mean μ₁(LC_k,pol_p) obtained from the distribution calculated from the range of temperatures lower than −3 °C, determines the threshold over which the surface is considered as being severely frozen. Thus, if ${\mathsf{\Delta }\mathsf{\sigma }}^0$ ≥ μ₁(LC_k,pol_p), the surface is assumed to be severely frozen (50% of the area of the normal function). The part where the normal function that corresponds to ${\mathsf{\Delta }\mathsf{\sigma }}^0$< μ₁(LC_k,pol_p) designates surfaces that are either unfrozen or having a mild-to-moderate freezing episode (e.g., plots sheltered from the cold). The mean value μ₂(LC_k,pol_p), obtained from the distribution using the temperature ranges [−3 °C, 0 °C], determines the threshold below which the surface is supposed to be unfrozen (i.e., ${\mathsf{\Delta }\mathsf{\sigma }}^0$ < μ₂(LC_k,pol_p)). Finally, the zone that corresponds to μ₂(LC_k,pol_p) ≤ ${\mathsf{\Delta }\mathsf{\sigma }}^0$ < μ₁(LC_k,pol_p) is considered as having a mild-to-moderate freezing episode. Figure 4 shows the flowchart of this algorithm.

3.4. Sentinel SAR Data Preparation

For both study sites and each date, Sentinel-1 SAR normalized backscatter values (${\sigma }^0$) were calculated at plot scale in both VV and VH polarizations by averaging ${\sigma }^0$ of all the pixels inside each plot, and for each polarization separately. For each polarization, the Sentinel-1 acquisitions from the morning passes (descending mode) and evening passes (ascending mode) were analyzed separately due to diurnal variations between both acquisitions. The diurnal variation is a result of the difference in the vegetation water content (VWC) between the morning and the evening. This difference in VWC could cause high difference in the radar-backscattering signal over vegetated plots between the morning and the evening acquisitions. In fact, several studies have reported that ${\mathsf{\sigma }}^0$ in the morning overpass registers higher values than ${\sigma }^0$ in the evening overpass [61,62]. Therefore, it is advisable to investigate separately each SAR temporal series acquired in the morning and the evening.

4. Results

4.1. Calculation of Thresholds for Frozen Detection

As described in Section 3.3, the thresholds over $\mathsf{\Delta }{\sigma }^0$ for the detection of frozen surfaces were obtained separately for each land cover class (three LCs) and each S1 polarizations (VH and VV). Figure 5 and Figure 6 show the distributions of $\mathsf{\Delta }{\sigma }^0$ that allowed the calculation of the thresholds after fitting normal functions. The distributions of $\mathsf{\Delta }{\sigma }^0$ were plotted for each LC using Sentinel-1 SAR images at two temperature ranges: −3 °C ≤ t ≤ 0 °C and t ≤ −3 °C. Then, from each fitted normal function, the mean was calculated (

Table 1. Calculated thresholds for each land cover class (LC) in VV and VH polarizations.

Pol.	Unfrozen $\mathsf{\Delta }{\mathit{\sigma }}^0\left(\mathsf{dB}\right)$			Mild-to-Moderate Freezing $\mathsf{\Delta }{\mathit{\sigma }}^0\left(\mathsf{dB}\right)$			Severe Freezing $\mathsf{\Delta }{\mathit{\sigma }}^0\left(\mathsf{dB}\right)$
	LC1	LC2	LC3	LC1	LC2	LC3	LC1	LC2	LC3
VH	$\mathsf{\Delta }{\sigma }^0$ < 3.5	$\mathsf{\Delta }{\sigma }^0$ < 2.8	$\mathsf{\Delta }{\sigma }^0$ < 2.1	3.5 ≤ $\mathsf{\Delta }{\sigma }^0$ < 5.3	2.8 ≤ $\mathsf{\Delta }{\sigma }^0$ < 3.5	2.1 ≤ $\mathsf{\Delta }{\sigma }^0$ < 2.9	$\mathsf{\Delta }{\sigma }^0$ ≥ 5.3	$\mathsf{\Delta }{\sigma }^0$ ≥ 3.5	$\mathsf{\Delta }{\sigma }^0$ ≥ 2.9
VV	$\mathsf{\Delta }{\sigma }^0$ < 2.5	$\mathsf{\Delta }{\sigma }^0$ < 1.7	$\mathsf{\Delta }{\sigma }^0$ < 1.6	2.5 ≤ $\mathsf{\Delta }{\sigma }^0$ < 4.0	1.7 ≤ $\mathsf{\Delta }{\sigma }^0$ < 2.2	1.6 ≤ $\mathsf{\Delta }{\sigma }^0$ < 2.4	$\mathsf{\Delta }{\sigma }^0$ ≥ 4.0	$\mathsf{\Delta }{\sigma }^0$ ≥ 2.2	$\mathsf{\Delta }{\sigma }^0$ ≥ 2.4

). The results show that the threshold values vary greatly depending on the land cover class and the polarization used. Table 1 shows that for LC1, which corresponds to cereals, severe freezing occurs when the calculated $\mathsf{\Delta }{\sigma }^0$ in VH polarization is greater than 5.3 dB, 3.5 dB for meadows (LC2), and 2.9 dB for orchards and vineyards (LC3). In addition, for VH, mild-to-moderate freezing occurs for $\mathsf{\Delta }{\sigma }^0$ between 3.5 and 5.3 dB for LC1, 2.8 and 3.5 for LC2, and between 2.1 and 2.9 for LC3. Finally, surface soils are considered as unfrozen if $\mathsf{\Delta }{\sigma }^0$ in VH polarization is lower than 3.5, 2.8 and 2.1 dB for LC1, LC2, and LC3, respectively. The thresholds in VV polarization, as shown in Table 1, are significantly lower in comparison to VH across the three land cover classes. For example, the threshold at which severe freezing occurs in VV is lower by 1.3 dB for both LC1 (4.0 vs 5.3 dB for VH), and LC2 (2.2 vs 3.5 dB for VH) in comparison to VH, and by 0.5 dB for LC3 (2.4 vs 2.9 dB for VH). The thresholds in VV polarization are also lower than those obtained for VH for the detection of mild-to-moderate freeze. For mild-to-moderate freeze, VV is lower than VH by roughly 1 dB for LC1 (2.5 vs 3.5 dB for VH) and LC2 (1.7 vs 2.8 dB for VH), and by 0.5 dB for LC3 (1.6 vs 2.1 dB for VH). These results are in accordance with the results obtained in the study of Baghdadi et al. [20], and confirm the higher potential of VH in comparison to VV for the detection of freezing episodes at the plot scale. It is worth mentioning that choosing lower threshold values, as is the case with VV polarization, would lead to a greater number of false detections of freezing episodes over each plot. However, after applying the temperature filter, some of these false detections could be removed.

4.2. Detection of Frozen/Unfrozen Plots

This section shows the results that correspond to only the first phase of the freezing detection algorithm over individual plots (i.e., without applying the temperature filter) in order to assess its performance. The algorithm has been tested on both study sites using all S1 images (acquired in Ascending and Descending modes, in both VH and VV polarizations). In our interpretation of the results, frozen plot detection is assumed plausible at a certain date, if the reported temperature by the nearest weather stations for each study site is below 3 °C. In reality, freezing only occurs when the temperatures are at, or below, 0 °C. However, it is possible that the nearest weather station to the studied site reports a higher temperature by a few degrees Celsius, while the temperature in a part of the study site is a few degrees Celsius below the recorded value (for reasons of distance or the morphology of the study site). Alternatively, false freeze detection corresponds to plots detected as frozen at a certain date while the reported temperature was above 3 °C, henceforth referred to anomalies. Finally, while we report the results of detections for both descending and ascending modes combined, the detection algorithm was run separately on each mode, and for each polarization.

4.2.1. Freezing Detection Results over Site A

Figure 7 shows the detection of frozen episodes over three plots (one for each land cover class) in site A (near Paris, France), using VH and VV polarizations. For the three land cover classes, we observe that the majority of detected freezing episodes across all land cover classes happened over the period that extends from November to February, which corresponds to the coldest months of the year, with temperatures (t) sometimes reaching 0 °C or lower. There was also one freezing episode that happened mid-April where the temperature dipped below 0 °C (t = −1.2 °C), and another potential freezing episode in late March with a temperature close to 0 °C (t = 2.9 °C). The majority of the anomalies were detected either before the freezing period (October), or at its end (mid-February onwards) and this was observed for both VH, and VV polarizations. These anomalies are caused by a combination of rainfall events occurring days or even weeks before the detection date of the anomalies, and the vegetation growth period, which starts in February. Irrigation events can equally contribute to the detection of anomalies; however, their effect remains much weaker since irrigation events affect some plots only from April. In addition, Figure 7 shows that the anomalies are more frequent in VV that in VH.

Over the entire plots in Site A, the percentages of frozen plots for each Sentinel-1 acquisition date and each land cover class (LC) are presented in Figure 8 for VH and VV polarizations. Figure 8a–c show that the true freezing episodes in VH polarization occurred between the period extending from November to mid-February, and then from late March until mid-April, and for the three land cover classes.

The results also show that the detected freezing episodes are strongly correlated with the temperature. Indeed, five of a total of eight dates that were detected as having a freezing episode had temperatures below 0 °C. Of these five freezing episodes, there are two dates (mid-November and mid-December) where the percentage of frozen plots was below 50% across the three land cover classes. For example, LC1 and LC2 report less than 20% of frozen plots for mid-November. For temperatures between 0 and 3 °C, only three freezing episodes were detected, with one occurring in early November (t = 0.4 °C) where the percentage of freezing plots was higher than 75% for LC1 and LC2, and 39% for LC3. For the two remaining freezing episodes, which occurred in late February and late March, the percentage of freezing plots was below 40% for LC1, and LC2, and below 60% for LC3, and the temperature averaged 2.9 °C. It is worth noting that while the results of LC3 are somewhat inconsistent with the findings for LC1 and LC2, this could be due to the low number of plots corresponding to LC3 (29 plots against 793 plots for LC1 and 179 plots for LC2) in site A, and also due to the lower thresholds used for LC3 in comparison to LC1, and LC2 (Table 1). During the freezing period ascertained through low temperatures (November through February), the anomalies are less frequent. For LC1, LC2, and LC3 in VH polarization (Figure 8a–c), there was only one date (27 January) which had, respectively, ~37%, ~17%, and 36% of plots falsely detected as frozen. There are four other dates with false detections in LC3 (three in November and one around mid-December); however, the percentage of plots with anomalies is less than 16%.

The percentage of plots detected as having anomalies at periods when there was freezing occurrence (temperatures higher than 3 °C) is in most cases lower than 50% for the three land cover classes (October, and March onwards). As mentioned before, the increase in the number of anomalies in the time period from March onwards is due to rainfall or irrigation events, while the anomalies detected in October are due to the moderate rainfall which occurred on 23 September 2018 (Figure 8, graph at the top). For this same period, part of the anomalies could be due to tillage, because agricultural plots, particularly class LC1, are plowed, which increases the radar signal and sometimes leads to an anomaly (especially when plowing leads to high soil roughness). These three main causes of anomalies (rainfall, irrigation and roughness), as explained in Section 3, cause some unfrozen surfaces to be falsely detected due to the high increase in the SAR signal in the S1 time series at irrigation/rainfall/tillage time, and therefore, we will have a significant Δσ⁰ since the σ⁰ ($d_n$), which corresponds to a $d_n$ close to a date of rain, irrigation or tillage, will most probably have a high value. However, as we will see in the upcoming sections, after the application of the temperature filter, these anomalies will be discarded.

Freezing episodes in the VH ascending mode (afternoon) for the three land cover classes (Figure 8a–c), are less frequent, due to the higher temperatures during the day. Therefore, only two freezing episodes where detected for LC1, and LC2 on 23 and 24 January (temperature ≈ 0 °C). For LC1, the percentage of frozen plots at these two dates is higher than 93% and ~60% for LC2. LC3 detected five freezing episodes, two of which occurred on 23 and 24 January, with the percentage of frozen plots being, respectively, 77.5% and 89.1%; the other freezing episodes were detected on 12 December (t = −2.3 °C), with less than 60% of plots detected as frozen, and two in late January/ early February 2019 (t ~ 2.65° for both dates), with less than 31% of plots detected as frozen. The anomalies detected using the SAR ascending VH acquisitions for LC1 and LC2 mostly follow the same dates as those found in the morning (i.e., minimal detection of anomalies during the freezing period; Figure 8a,b) and the majority of detected anomalies outside of this period do not surpass 50% of the number of plots. The high number of detected anomalies in LC3 (Figure 8c), as we previously mentioned, is not indicative of the bad detection of our algorithm for LC3 due to the low number of plots in site A, and also due to the lower thresholds used for LC3, in comparison to LC1, and LC2 (Table 1).

The freezing detection using SAR acquisitions in VV polarization often shows similar performance in comparison to VH polarization for the morning period (Figure 8d–f). For example, on 3 November and 12 February, where the reported temperatures were, respectively, ~0 and ~−1 °C, our algorithm applied over S1 images produced similar results as VH for the three land cover classes. However, VV was able to detect slightly more freezing plots on 27 March in comparison to VH (more than 70% for VH and less than 40% for VH); however, this is unlikely to be true since the reported temperature at these dates is around 3 °C. The higher rate of freezing detection at these dates using VV in comparison to VH is due to the lower threshold values used for VV in comparison to VH (Table 1), which causes the algorithm to detect more freezing episodes with lower thresholds. Indeed, this is also noticeable with the detection of freezing episodes in the afternoon, where, using VV, the algorithm was able to detect slightly more plots in comparison to VH. The detection of anomalies using VV polarization follows the same pattern as VH polarization. Indeed, the algorithm detected minimal anomalies during the period between November and February, and most of the detected anomalies appear after rainfall events. However, also due to the lower threshold used for VV, our algorithm detected a higher percentage of plots with anomalies in comparison to VH.

4.2.2. Freezing Detection Results over Site B

Figure 9 shows the performance of our algorithm across three reference plots corresponding to three land cover classes over site B (near Strasbourg, France), using S1 images in VH (Figure 9a–c) and VV (Figure 9d–f) polarizations.

Over the reference plot corresponding to LC1 (Figure 9a), 13 freezing episodes (morning and afternoon combined) between mid-November and mid-February were detected using VH polarization. In the same time period, the algorithm did not detect any anomalies. For the reference plot corresponding to LC2 (Figure 9b), 19 freezing episodes (morning and afternoon) were detected in VH polarization also between mid-November and mid-March. Of these 19 freezing episodes, two episodes in early February corresponded to anomalies because the temperature at the time of acquisition of S1 was above 5 °C for both dates, which suggests a very low to no probability of freezing. Finally, for LC3 (Figure 9c), the algorithm detected 13 freezing episodes (morning and afternoon) for the period between mid-November and mid-February, with only one anomaly in mid-November. The majority of the anomalies for the three land cover classes appear to be detected outside of the freezing period that extends from November to February; this is in accordance with the detected anomalies in site A. Moreover, the higher number of detected freezing episodes over site B in comparison to site A is mostly due to the lower yearly average temperatures in comparison to site A, which increases the occurrence of freezing events.

Our algorithm detected fewer freezing episodes across the three land cover classes when used with VV polarized S1 images over the same reference plots. For LC1 (Figure 9d), only six freezing episodes were detected for the period between mid-January and mid-February, and one freezing episode was detected in early March. However, the algorithm also detected seven anomalies between mid-February and mid-March. For LC2 (Figure 9e), the algorithm detected ten freezing episodes between mid-November and mid-February in VV polarization with no anomalies. Finally, for LC3, the algorithm detected 12 freezing episodes between mid-November and mid-February, of which two were anomalies (temperature above 5 °C for both dates).

Figure 10 shows the distribution of the frozen episodes and anomalies, across all the plots in each land cover class for both VH and VV acquisitions in ascending and descending modes. The period of the freezing episodes’ occurrence is similar to those shown for site A, with the overwhelming majority of freezing episodes detected between mid-November and mid-February, and this was observed for both VH and VV polarization (Figure 10).

The percentage of frozen plots is tightly coupled with the temperature for both VH and VV polarizations. For temperatures below 0 °C, and using VH polarization, the percentage of detected frozen plots for LC1 was higher than 50%, and 80% of plots were detected as frozen in four out of the five freezing episodes. For temperatures around 0 °C (t between 0 and 1 °C), at least 50% of plots in LC1 were detected as frozen in 50% percent of freezing episodes. Finally, for temperatures above 1 °C, the percentage of frozen plots was less than 15%. For LC2, and LC3, in VH polarization, similar patterns to LC1 were observed. However, for the lowest temperatures (t < −1 °C), our algorithm detected fewer frozen plots in LC2, and LC3 in comparison to LC1 (Figure 10c,d). For the ascending (afternoon) acquisitions in the VH polarization, the highest percentage of frozen plots were observed on 21 January and 24 January 2019 for all the land cover classes, with a reported temperature of 1.5 °C. On 21 January 2019, the percentage of frozen plots was 91% for both LC1 and LC2, and 94.5% for LC3. Alternatively, on 24 January 2019, the percentage of frozen plots was less than the previous date, with, respectively, 75%, 32.5%, and 51.4% for LC1, LC2, and LC3.

Using VH polarization, the detected anomalies for LC1, LC2, and LC3 were very few during the freezing period (November to mid-February) and for both descending (morning), and ascending (afternoon) acquisitions. In addition, the anomalies during this period, when detected, did not exceed 5% of the total number of plots for the three classes. However, as with site A, the detected anomalies mostly occurred in October due to three heavy rainfall episodes happening in September, and starting from mid-February onwards, and these are due to a combination of rainfall, irrigation, tillage and/or vegetation growth.

The freezing detection results across the three land cover classes using VV polarization (Figure 10d–f) show mostly the same detection capabilities as those obtained using VH polarization for the period between November and mid-February, and for both descending (morning), and ascending (afternoon) acquisitions. However, due to the lower threshold values used for VV, our algorithm detected more frozen plots on 13 March 2019 in comparison to VH, and also a lot more anomalies starting from mid-February onwards.

4.3. Improved Freezing Detection by Applying the Temperature Filter

The results obtained in the previous section show that our algorithm was successful in detecting frozen plots with a minimum number of anomalies during the freezing period between November and February. These results are consistent across both study sites, and in both descending (morning), and ascending (afternoon) acquisition modes. However, for the period starting from mid-February, and for the month of November, there were many more anomalies detected (detections with temperature > 3 °C). To remove these anomalies, the second phase of our algorithm consists of applying a filter that removes all detections when the reported temperature from the nearest weather stations is above 3 °C. This section also presents the distribution between mild-to-moderate freezing and severe freezing (depending on threshold values as shown in Table 1).

Figure 11 shows the distribution of the severity of freezing in site A across the land cover LC1, LC2, and LC3 in both VH (Figure 11a–c) and VV (Figure 11d–f), respectively. The results indicate that the severe freezing episodes occurred exclusively when the temperatures were below 0 °C for both VH and VV polarizations in descending mode (morning). For the VH polarized images, the percentage of severe frozen plots was at its highest when the temperature was at its lowest (~−2.5 °C on 26 December 2018), with the percentage of severely frozen plots roughly equal to 40%, 77 % and 75% for LC1, LC2, and LC3, respectively. The second date with the most severe frozen percentages occurred on 14 April 2019 (temperature ~−1.5 °C) with a percentage of severe freezing of ~45% for both LC1 and LC2, and ~62% for LC3. In addition, the detection results using VV polarized images also showed that the highest percentage of severely frozen plots also occurred on 26 December 2018, and 14 April 2019. Finally, our algorithm using both VH and VV polarized images did not detect any severely frozen plots in the afternoon period.

The results for site B (Figure 12) show that the detection of the severe freezing also occurred on the dates with the lowest temperature, and this is in accordance with the results obtained in site A. For VH and VV polarized images in the descending mode (morning), the dates with the most detected severe frozen plots occurred when the reported temperatures were at their lowest (t < −3 °C on 23 and 24 December 2018), and for the three land cover classes. On 23 December 2018, the percentage of severely frozen plots was ~60% for LC1 in both VH and VV, ~65% and ~68% for LC2 in, respectively, VH and VV, and ~50% for LC3 in both VV and VH. The percentage of severe frozen plots was slightly higher on 24 December, with ~70% of plots detected as severely frozen for LC1 in both VH and VV, ~68% and ~76% for LC2 in, respectively, VH and VV, and ~76% for LC3 in both VH and VV. Finally, for the images acquired in both VH and VV in the ascending mode (afternoon), our algorithm did not show any occurrence of severe freezing across the three land cover classes.

4.4. Comparison with ERA5-Land Dataset

The assessment of accuracy for our algorithm to detect frozen/unfrozen plot surface states is a challenging issue, as direct ground measurements of the actual plot surface states are not available. Therefore, it is necessary to resort to proxy variables, such as the air temperature [60] or soil temperature [28]. However, proxy variables might not necessarily represent the surface state and may sometimes introduce uncertainties [16]. In addition, the proxy variables (soil and air temperatures) should be acquired in the same timeframe as Sentinel-1 acquisitions, as any interpolation of temperature data to Sentinel-1 acquisition time would inevitably lead to uncertainties. In this study, hourly soil and air temperature data from ERA5-Land were used. ERA data products have been utilized in a multitude of studies for the validation of freeze/thaw products [15,16]. A previous ERA product, namely the ERA Interim, was found to be the most accurate reanalysis data product available [49,50]. Moreover, ERA5-Land, which reports more accurate air temperatures than ERA interim [63], has a standard deviation to in situ air temperatures of 1.4 °C [64]. Over both of our study sites, the reported standard deviation to the in situ air temperatures was ~1 °C. Finally, due to the large difference in spatial resolutions between ERA5-Land data products (9 Km) and our results (plot scale), a qualitative analysis will therefore be made between ERA5-Land (soil and air temperatures) and the percentage of frozen/unfrozen plots at a certain date, using the three land cover classes combined (i.e., percentage of the combined detected mild-to-moderate and severe frozen plots over the entire study site).

Over site A (near Paris), Figure 13a shows that when both the soil and air temperatures reported by ERA5-Land were below 0 °C, the highest percentage of frozen plots detected using VH occurred for both the both morning and afternoon acquisitions. For these temperatures, the reported percentage of frozen plots was between 42% and 93%, with 85% of frozen episodes having more than 75% of plots detected as frozen. For ERA5-Land temperatures (both air and soil temperatures) higher than zero (between 0 and 3 °C, since the temperature filter is applied for t > 3 °C), the reported percentage using VH was lower, and varied between a few plots detected as frozen and approximately 40%. Freezing detection results obtained from VV acquisitions (Figure 13b), as seen in the previous sections, show the same behavior in comparison to VH.

The comparison of our results to ERA5-Land’s soil and air temperatures for site B (Strasbourg) show that the detection results obtained using VH and VV are mostly similar (Figure 14a,b). For both polarizations, as in site A, the most detected frozen plots occurred on dates when both ERA5-Land air and soil temperatures were significantly lower than 0 °C (between 19 and 25 January 2019 and 4 and 12 February 2019). On these dates, the percentage of frozen plots varied between 62% and 95%. For less cold ERA5-Land soil temperatures (0–1 °C), the percentage of frozen plots decreased; with a percentage of detected frozen plots ranging from 15% to 67% with a mean of 33% using VH, and from 11% to 53% with a mean of 59% using VV. For ERA5-Land soil temperatures between 1 and 3 °C, the percentage of frozen plots did not exceed 7.4% using VH with a mean of 2.5%. The freezing detection results using VV were similar to VH for ERA5-Land soil temperatures between 1 and 3 °C. However, on 12 March 2019 (morning), when both ERA5-Land soil and air temperatures were roughly 1.3 °C, our algorithm detected 29% of plots being mild-to-moderately frozen using VV against ~7% using VH. Finally, on 31 January 2019, the reported soil temperature by ERA5-Land at Sentinel-1 acquisition time (~6 PM) was ~0 °C, but our algorithm did not detect any plots as frozen using either VH or VV. However, on this date, the reported ERA5-Land air temperatures from several hours before S1 acquisition, and up-to S1 acquisition time (~6 P.M.), was relatively high (~3 °C). Therefore, we suspect, given the uncertainties on ERA5-Land soil temperatures and the warm air temperatures, that no freezing occurred on this date and time.

5. Discussion

5.1. Freezing Detection Using Sentinel-1 σ° SAR Backscattering

As evidenced in Section 4.1, Section 4.2 and Section 4.3, the behavior of the S1 SAR backscatter response differed greatly according to the land cover class. Therefore, it was necessary to detect freezing episodes with different thresholds based on the investigated land cover classes and the S1 polarization. In this study, the proposed approach to calculate the threshold values for freezing detection using S1 SAR images in both polarization modes (VH and VV) appeared to perform relatively well in determining the best thresholds for each land cover class.

The detection of freezing events was performed by calculating the difference between the investigated ${\sigma }^0$ at date ‘$d_i$’, and a mean value of the S-1 SAR signal obtained by using the value of three previously detected maximum ${\sigma }^0\left(max\_{\sigma }^0\right)$ from several prior days. $max\_{\sigma }^0$ are obtained by detecting the maximum S1 SAR backscattered coefficients in a moving window of 15 days from and up-to date ‘$d_i$’. We chose the mean of the last three $max\_{\sigma }^0\left(mean\_{\sigma }^0\right)$, instead of just the latest one, for several reasons. In this study, we rely on a change detection approach (difference between the acquired ${\sigma }^0$ at date ‘$d_i$’ and $mean\_{\sigma }^0$) to determine a freezing occurrence. Thus, using only the latest $max\_{\sigma }^0$ may make our model susceptible to making false detections due to several phenomena (vegetation growth, rainfall, irrigation, tillage, or even signal noise) that lead to an increase in the difference between $max\_{\sigma }^0$ and ${\sigma }^0$ at date ‘$d_i$’. In addition, such phenomena can also lead to uncertainties when distinguishing between mild-to-moderate and severe freezing episodes. For example, after a rainfall or irrigation event, and even in less cold temperatures (between 0° and 3 °C), relying on only the latest $max\_{\sigma }^0$, our algorithm can sometimes falsely detect mild-to-moderate frozen episodes as severe. Therefore, to mitigate the effect of these phenomena (vegetation growth, irrigation, rainfall, and tillage), the use of $mean\_{\sigma }^0$ would improve our prediction accuracy. Indeed, our algorithm did not detect any presence of severe frozen plots with less cold temperatures in both the morning and afternoon periods. Moreover, severely frozen plots were detected only when the temperatures were at their lowest (period between November and February), and this was observed over both study sites. For this period (between November and February), the number of falsely detected frozen plots was also greatly reduced.

Finally, the results are in accordance with the findings obtained via simulation in the study of Baghdadi et al. [20]. In their study, they also showed a strong potential for Sentinel-1 data to map freezing episodes at the plot scale in an operational capacity. The thresholds obtained are also consistent with the findings of the simulations. Simulation results in Baghdadi et al. [20] found that a signal drop of at least 3-4 dB corresponded to a freezing episode (depending on the type soil). With actual Sentinel-1 data, we considered for LC1 that freezing occurs when the signal drops below 3.5 dB, 2.8 dB for LC2, and 2.1 dB for LC3, and a severe freezing is detected for signal drops higher than 5.3 3.5 and 2.9 for, respectively, LC1, LC2, and LC3.

5.2. Effect of the Temperature Filter

While false detections during the freezing period over our study sites (between November and February) did not exceed most of the time 10% across the land cover classes; false detections were mostly present outside of this period, and especially after rainfall events (November, and February through April), and a combination of rainfall, and irrigation events from May onwards. Therefore, the inclusion of in situ air temperature filter was found to be particularly relevant to entirely discard these false detections. In this study, we chose an air temperature threshold of 3 °C, over which we removed any freezing detection. In reality, freezing occurs when soil temperature, not air temperature, is equal to or below 0 °C. However, due to several reasons—(1) the tight coupling between air temperatures and that of top soil [65]; (2) The distance between our study sites and their corresponding weather stations and the difference in morphologies; (3) the fact that soil temperatures are generally lower than air temperatures [17]—we believe an air temperature threshold of 3 °C would have a better likelihood of the air temperature representing actual frozen conditions [60].

5.3. Freezing Detection Approach

Overall, the freezing detection results of our algorithm using either VH or VV were very similar over both sites, and across the three land cover classes. In essence, there were no freezing episodes that were detected using a particular polarization that were not detected in the other. Moreover, when we detected a freezing episode on a particular date, we were able to detect it for the three land cover classes. However, we suspect that prediction results using VH should present better accuracies, as the threshold ranges are higher in VH than VV. While both VH and VV falsely detected freezing episodes, mostly on the same dates, the number of falsely detected frozen plots was less prominent in VH than in VV, for both morning and afternoon acquisitions, and in both study sites. Moreover, the higher threshold values calculated for LC1 in comparison to LC2, and LC3 showed that LC1 is, on occasion, more precise when detecting frozen plots over both study sites. Indeed, for a given freezing episode, the number of detected plots in LC1 was the highest, followed by LC2 and LC3.

5.4. Comparison to ERA5-Land Soil and Air Temperatures

The qualitative analysis between the detection results using our algorithm and ERA5-Land’s air and soil temperatures supported our previous findings. In effect, our algorithm is capable of detecting the same freezing episodes detected using ERA5-Land. However, while both product performances (our algorithm and ERA5-Land) were similar, the qualitative precision of our frozen maps appears to be higher. Using our algorithm with the high spatial resolution S1 images allowed us to better differentiate the local variations in freezing episodes at the plot scale in comparison to ERA5-Land, which is unable to detect them due to its low spatial resolution (9 km). For example, on 27 December 2018 and 25 February 2019, in site A, the reported soil temperature by ERA5-Land was, respectively, −2.7 and 2.6 °C, which suggests that two freezing events with different intensities occurred. However, Figure 15a,b, show that our product is more capable of describing the intensities of these freezing episodes and discerning which plots were affected. Similarly, over site B, on 22 January 2018 (t = −5.5 °C) and 25 January 2019 (t = 0.9 °C), we observe that for the lowest temperatures, our product detected most of the plots as frozen, with almost half detected as being severely frozen (Figure 15c). Alternatively, on 22 January 2019, with a reported temperature of 0.9 °C, our algorithm did not detect any severely frozen plots, and only slightly more than half of the plots were detected as being mild-to-moderately frozen (Figure 15d).

5.5. Frozen Mapping Qualitative Analysis

Even if Figure 15 shows that the global mapping of frozen surfaces is coherent with the temperatures, some plots are detected as unfrozen with surrounding severely frozen plots. On 27 December 2018, 19% of plots in Site A (Figure 15a) were detected as unfrozen, while the temperatures showed favorable freezing and severe freezing conditions. On this date, 43% of the plots were detected as being moderately frozen, and 38% as severely freezing. While it is not uncommon for some plots to actually be unfrozen at a temperature of −2.7°C, the question pertaining to the relevance of mapping unfrozen plots around severely frozen plots arises. For all the plots detected as unfrozen, the radar signal on December 27 has decreased significantly but did not cross the calculated thresholds in order to be detected as frozen. In fact, 50% of the plots detected as unfrozen in LC1 registered a drop in the radar signal between 2.8 and 3.48 dB (at most 0.7 dB from the chosen freezing threshold, which was set at 3.48 dB). Moreover, 34% of the plots detected as unfrozen in LC1 registered a drop in the radar signal between 2.3 and 2.8 dB (at most 1.2 dB from the selected threshold). As for LC1, the same patterns also appeared for LC2, and LC3. For LC2, and LC3, 50% of the plots detected as unfrozen were only 0.7-0.8 dB below the chosen detection threshold values, and 34% were 1.3–1.4 dB below the thresholds. Therefore, only 16% of the plots mapped as unfrozen on 17 December 2018 have a weak signal decrease. These plots, which correspond to only 3% of the total number of plots within site A, could be falsely detected as unfrozen, most probably due to variations in soil roughness (plots tilled a few days before the mapping date, which causes an increase in soil roughness and consequently of the radar signal). Other factors, such as the soil type or the presence of hedges, can also explain in part these behaviors. The vegetation effects are less probable to cause any misdetections due to the low vegetation growth at this time period.

On 22 January 2019, 12% of the plots within site B (Figure 15c) were detected as being unfrozen, while the temperatures were well below −5 °C, which suggests severe freezing conditions. On this date, 26% of the plots were detected as being moderately frozen and 62% as being severely frozen. Overall, 50% of the plots in LC1 which were detected as unfrozen showed a drop in the radar signal between 2.5 and 3.48 dB. However, since the calculated threshold for LC1 was 3.48 dB, these plots could not be classified as frozen (the radar signal drop was at most 1 dB below the threshold). Moreover, 34% of the plots detected as unfrozen in class LC1 registered a decrease in the radar signal between 1.7 and 2.5 dB (at most 1.8 dB below the selected threshold). For LC2 and LC3, 50% of the plots detected as unfrozen were 0.7 dB below the calculated thresholds, and 34% were below the calculated thresholds by 1.2 to 1.3 dB. As such, only 50% of the plots which were detected as unfrozen on 22 January 2019 had a radar signal decrease below the detection threshold of less than 1 dB for LC1, and of less than 0.7 dB for LC2 and LC3, which corresponds to around 6% of the total number of plots within site B. These weak drops in the radar signal for certain plots could be explained, as in site A, by an increase in the soil roughness several days before the acquisition of 22 January 2019. Finally, the results presented in Figure 15b,d, do not show any aberrant detections, as we only detected unfrozen to mild-to-moderately frozen plots, which is consistent with temperatures around 0 °C.

5.6. Strengths and Limitations

The results of our algorithm show that freezing episodes could be successfully detected using the high-resolution Sentinel-1 SAR images, with, at worst, a temporal resolution of six days. However, the implementation of the proposed method in a near real-time scenario primarly depends on the delivery time of S1 images. S1 images are usually delivered by ESA in a “fast 24 h” delivery mode. This mode ensures that S1 images are available for download 24 h after the satellite acquisition. Considering that the pre-processing of S1 images and applying the proposed algorithm is automatically performed with minimum human intervention, frozen plots could be detected in a few hours after the reception of S1 images based on the the size of the studied area.

To detect a freezing episode, we rely on a simple but effective algorithm that requires only the calculation of two thresholds (to distinguish between unfrozen, mild-to-moderately frozen, and severely frozen) for each land cover class, and two additional datasets (land cover classes, and in situ air temperatures). The land cover classes were provided from the detailed graphical parcel registry (RPG), which is available over all France. However, RPG was only needed to derive three super land cover classes (LC), where each LC encompassed several crop types. Hence, if no land cover maps are available, the user can calculate a single set of thresholds for the entirety of the studied plots, which might lead to uncertainties in the detection of freezing in some of the land cover classes. In addition to the land cover classes, we used in situ air temperatures from the nearest weather station to (1) calculate the threshold-values, and (2) apply the temperature filter. Therefore, in remote areas, or when no in situ temperature data are available, they can be substituted by the less accurate, but globally available, ERA5-Land’s air or soil temperatures. Using these three datasets (Sentinel-1 SAR, land cover classes, air temperatures), our algorithm provided near real-time results of plots’ frozen states, which distinguishes between unfrozen, mild-to-moderately frozen, and severely frozen.

The strength of our algorithm is multi-faceted. First, it is computationally inexpensive, and has the potential to detect freezing episodes at near real-time. In general, S1 SAR images are available 24 h after acquisition; therefore, frozen plots could be detected a day after the event. Secondly, over agricultural areas, vegetation growth, tillage, and irrigation can cause large fluctuations in the backscatter response, which has the adverse effect of causing uncertainties when detecting freezing episodes, as witnessed in previous studies [15,16]. In this study, we do not rely on the direct value of the ${\sigma }^0$ to detected freezing, but rather on the difference between the acquired ${\sigma }^0$ and the mean of the previous three calculated $max\_{\sigma }^0$. This approach mitigated the effects mentioned previously, and allowed us to precisely detect freezing episodes over agricultural areas and pastures.

There are, however, limitations in the operational large-scale mapping of freeze/thaw detection at the plot scale. First, the validation of freeze/thaw maps is strewn with pitfalls given the scarce availability of relevant data such as soil ice contents. High-resolution meteorological data are also not free, or available as open access, and the only accessible temperatures data, as mentioned in our article, are several km away from the study sites. Therefore, in this study, as many previous studies (e.g., [15,16,17,28,47]) that attempted to detect freeze/thaw episodes, we relied on proxy variables, such as air and soil temperatures from ERA5-Land, in order to assess the performance of our algorithm. For example, Naeimi et al. [16] used ERA interim temperatures in order to derive a freeze/thaw product and thus a quantitative analysis was possible. However, such quantitative analysis was not relevant for our study due to the following reasons:

• While ERA5-Land soil temperatures can give a good representation of a soil’s state, temperature is not the only parameter that causes a soil surface to freeze/thaw, and this is especially true for temperatures around 0 °C, during the afternoon period or during spring. Therefore, deriving a freeze/thaw map from ERA5-Land soil temperatures is not straightforward. In general, for temperatures lower than 0 °C, there are no issues when comparing our results and the reported temperatures by ERA5-Land. As seen in the results, for low temperatures (≤ 0 °C), the number of detected frozen plots was the highest. In addition, for slightly less cold temperatures (> 0°C), and while the number of detected plots was low, we have no way of determining if a freezing event occurred, and if so, the extent of its severity.
• The large difference between ERA5-Land’s 9km spatial resolution and our plots’ surface area (several hundred square meters) makes it very hard to do any meaningful statistical analysis. Even if we derive a freeze/thaw map using ERA5-Land soil temperatures using, for example, a threshold of 0 °C, below which a surface state is considered frozen, it is not possible to aggregate freeze/thaw detections from plot resolution to ERA5-Land’s resolution.

Given the presented uncertainties, in this study we opted to perform qualitative analysis instead of quantitative analysis. In fact, given the unavailability of free and open access of high-resolution temperature data, it is difficult to do otherwise. However, the use of freely available low-resolution in situ temperature data allows us to develop an operational methodology for the mapping of freezing episodes that is easily reproducible elsewhere, with the minimum amount of data.

Secondly, the proposed method is capable of showing plots that are supposed to be moderately to severely frozen. What is certain, even if we have no ground truth, and this is thanks to the simulation data present in the paper of Baghdadi et al. [20], is that if the signal drops for bare soils by at least 3 dB, frost has occurred. In contrast, the distinction between moderate and severe freezing remains arbitrary due to the chosen threshold that distinguishes between mild-to-moderate and severe freezing conditions. As a result, a small portion of the plots estimated to be severely frozen may actually be moderately frozen and a small portion of the plots mapped as moderately frozen may actually be severely frozen. What is certain, however, is the presence of multiple freezing conditions that we are capable of detecting. In any case, due to the lack of very high spatial resolution field data and using a threshold-based method, there may be a slight mixing between classes, and it is difficult to do otherwise for large-scale mapping.

A third limitation is linked to the temporal resolution of Sentinel-1 images. Sentinel-1 does not provide daily images, and thus frozen episodes that occur on days when there are no acquisitions cannot be detected. Sentinel-1 revisit times depend on the region, and in worse case acquisitions can occur once every six days.

Finally, the detection of freezing episodes in areas encountering seasonal freezing (frost persists for several months) could be problematic. In this study, we calculate the mean of previous $max\_{\sigma }^0$ by making the assumption that $max\_{\sigma }^0$ did not occur on a freezing date. However, in areas affected by seasonal freezing, the mean of the previous $max\_{\sigma }^0$ will not be guaranteed to be detected on an unfrozen date. Therefore, the difference between ${\sigma }^0$ and the mean of previous $max\_{\sigma }^0$ will be very small, and a false detection will be made.

6. Conclusions

In this study, we propose a new approach for mapping freezing episodes over agricultural areas, at plot scale, and in near real-time using Sentinel-1 SAR images acquired in VH and VV polarizations. Our approach relies on change detection in the S1 SAR backscattering coefficients (${\sigma }^0$) to detect frozen plots and can be divided into two major parts. The first part is responsible for calculating thresholds in order to distinguish between unfrozen, mild-to-moderately frozen, and severely frozen plots. Given that SAR backscatter fluctuations due to freezing can greatly vary between land cover classes, we opted to calculate the thresholds for each land cover class, and for each polarization (VH and VV). The second part of our algorithm is charged with the detection of frozen/unfrozen plots. The detection of a plot’s state at a particular date ‘$d_i$’ is done by calculating the difference between ${\sigma }^0$ at date ‘$d_i$’ and a mean value of the S-1 SAR signal obtained by using the value of three previously detected maximum ${\sigma }^0\left(max\_{\sigma }^0\right)$ from several days prior. $max\_{\sigma }^0$ values are obtained by detecting the maximum S1 SAR backscattered coefficients in a moving window of 15 days from and up to date ‘$d_i$’. Finally, due to several phenomena which can cause false freezing detection (rainfall, irrigation, tillage), a filter based on the registered air temperature at S-1 SAR’s acquisition time was introduced at the end of the detection chain. Our algorithm therefore will consider a particular plot as unfrozen if the reported temperature is higher than 3 °C.

Our approach was tested over two agricultural sites in France, site A (near Paris), and site B (near Strasbourg). Site A has, in general, higher mean yearly temperatures than site B, and the occurrence of freezing episodes is less common. The results show that our algorithm was capable of detecting all major freezing episodes in both sites, using VH or VV, and in both the morning and afternoon periods. The results also show that the detection of frozen plots using VH could be less ambiguous than using VV, given the higher thresholds used in VH than in VV (the radar signal decay due to freezing is greater in VH than in VV), which resulted in a lower number of false detections. The higher thresholds used for LC1, in comparison to LC2 and LC3, also showed that LC1 has, on occasion, better accuracy when detecting frozen plots over both study sites. Finally, the qualitative analysis between our results and ERA5-Land’s soil and air temperatures showed high agreement. In essence, there were no freezing episodes detected in ERA5-Land that were not detected using our algorithm. However, our algorithm, which uses Sentinel-1 high resolution SAR images, is capable of detecting freezing episodes at the plot scale, in comparison to ERA5-Land which shows a coarser representation of freezing episodes (spatial resolution of 9 km).