Potential of Sentinel-1 SAR to Assess Damage in Drought-Affected Temperate Deciduous Broadleaf Forests

Abstract

Understanding forest decline under drought pressure is receiving research attention due to the increasing frequency of large-scale heat waves and massive tree mortality events. However, since assessing mortality on the ground is challenging and costly, this study explores the capability of satellite-borne Copernicus Sentinel-1 (S-1) C-band radar data for monitoring drought-induced tree canopy damage. As droughts cause water deficits in trees and eventually lead to early foliage loss, the S-1 radiometric signal and polarimetric indices are tested regarding their sensitivities to these effects, exemplified in a deciduous broadleaf forest. Due to the scattered nature of mortality in the study site, we employed a temporal-only time series filtering scheme that provides very high spatial resolution (10 m $\times$10 m) for measuring at the scale of single trees. Finally, the anomaly between heavily damaged and non-damaged tree canopy samples (n = 146 per class) was used to quantify the level of damage. With a maximum anomaly of −0.50 dB ± 1.38 for S-1 Span (VV+VH), a significant decline in hydrostructural scattering (moisture and geometry of scatterers as seen by SAR) was found in the second year after drought onset. By contrast, S-1 polarimetric indices (cross-ratio, RVI, $H\alpha$) showed limited capability in detecting drought effects. From our time series evaluation, we infer that damaged canopies exhibit both lower leaf-on and leaf-off backscatters compared to unaffected canopies. We further introduce an NDVI/Span hysteresis showing a lagged signal anomaly of Span behind NDVI (by ca. one year). This time-lagged correlation implies that SAR is able to add complementary information to optical remote sensing data for detecting drought damage due to its sensitivity to physiological and hydraulic tree canopy damage. Our study lays out the promising potential of SAR remote sensing information for drought impact assessment in deciduous broadleaf forests.

1. Introduction

Droughts are becoming more severe and frequent as the Earth’s climate warms [1,2,3], exacerbating forest stress and tree mortality across all major types of forest ecosystems [4]. Accurate understanding of drought effects by continuously monitoring forest decline is, therefore, crucial given the central role of forest ecosystems for global biogeochemical fluxes, human welfare, and maintaining biodiversity and ecosystem services [4,5,6]. Assessing tree mortality in the field is challenging and costly, which sets limits to both spatial and temporal coverage of ground assessments. By contrast, remote sensing can provide large-scale and temporal high-resolution data suitable for forest damage assessments. However, understanding the role and dimension of water deficits and forest structural changes to tree mortality still needs to be better linked to satellite data.

In the optical domain of remote sensing, the relation to the reflectance of the leaf area [7], photosynthetic activity [8], and surface temperature [9] can be used to derive damage assessments and has been applied to a variety of sensors and frequency bands [10,11], including large-scale mapping [12]. Due to their short wavelengths, however, only the reflections of the leafy components of tree canopies were captured. On the other hand, since measurements at longer wavelengths have higher penetrations in the canopy, their signals are physically more related to phenomena inside the tree crown, consisting of leaves and woody compartments, each of which suffers a water deficit under drought conditions [13]. In the microwave domain, plants emit, absorb, and scatter energy, varying with the forest’s dielectric properties (their water content) and structures of plant compartments (distribution of leaves, branches, and stems), owing to the complex scattering mechanisms from forests [14]. By exploiting this relationship, microwave sensors hold great potential as damage-sensitive forest health monitoring instruments that can further provide spatially gapless and temporally dense observations due to their ability to penetrate clouds. However, the main challenge remains to quantify the sensitivity of the various microwave instruments to drought-induced damage to the canopy.

The estimation of water stress is a strong indicator of drought-induced damage. Good sensitivity for measuring water stress was recently found with microwave signals mainly for sensors with high temporal coverage, primarily space-borne scatterometers [15,16], radiometers [17,18,19], or tower-based experiments [20,21,22,23], e.g., due to exploiting symptoms of water stress in diurnal variations of water uptake. However, due to the coarse spatial resolutions of radiometers and scatterometers (typically > 10 km), the retrievals of plant hydraulic parameters are still limited to regional scales.

Tree damage and mortality are oftentimes underestimated when signals of affected and unaffected vegetation mix [6,24]. The damage patterns then remain outside the measurable monitoring scale when the sensor resolution is lower than the scale of single canopies. Here, synthetic aperture radar (SAR) can fill the observation gap. Operating at very high-resolution (10 m × 10 m), the ESA Copernicus twin Sentinel-1A/B (S-1) SAR satellites ($\lambda$ = 5.6 cm, f = 5.405 GHz) are vital tools for active microwave measurements. Due to their free and open data policies and dense observation time series with 12- (global) and 6-day (Europe) repeat cycles, the sensors provide global dual-polarimetric (VV-VH) C-band SAR data of high accessibility [25]. In comparison to longer L-bands (strong ground contributions) and shorter X-bands (mostly leaf contributions), C-band SARs penetrate substantial parts of the tree canopies while measuring the integrated water and structural properties [26]. The density of the S-1 time series is a game changer for the damage assessment of single canopies since now an effective temporal speckle filtering scheme (i.e., noise reduction) can be employed without losing spatial resolution [27].

In the field of forest ecology, S-1 data were used for forest type mapping [28], change detection in deforestation [29,30], windthrow mapping [31], retrieving phenological parameters [32,33], live fuel moisture content [34], and tracking insect defoliation [35,36]. Furthermore, S-1 could be used to derive vegetation water status in agricultural settings [37,38]. However, the sensitivity of temporally dense C-band SAR data still needs to be explicitly tested to assess signal sensitivity to drought-damaged deciduous broadleaf forest canopies.

As the microwave signal relates to water and structure, the backscatter coefficient is a mixed signal [39]. Given this complexity of absorption and scattering in tree canopies and the historical sparseness of SAR time series, damage assessments have long been ineffective at disentangling the plant water status, phenology/structure, and biomass variation [40]. At the timescales of weeks to months, and in the presence of strong droughts, the progressing dehydration of plant tissue becomes the primary driver of tree mortality [13]. The dense S-1 time series can here be used to filter/aggregate acquisition periods of relative constant biomass (growing season), smoothing over short-term influences (e.g., rain interception) so that signal changes can be related to water status variability (tissue dehydration), structural decay (e.g., early leaf loss) and, finally, drought-induced damage.

Steele-Dunne et al. [41] showed that soil, branch, and trunk influences on the C-band backscatter become negligible in broadleaf canopies when sufficiently wet, at least given the model presumptions on tree geometry of very dense aspen trees [42]. Generally, with increasing canopy water content, non-canopy C-band scattering becomes proportionally attenuated by the canopy volume. Still, in the event of droughts, leaves dry out and shed, thus leading to (1) a decrease of backscatter from leaves, (2) increasing scattering from branches, trunks, and soil, and (3) changes in scattering mechanisms from the multiple scattering of canopies to the single scattering of individual branches and twigs. S-1 polarimetric properties are used to test these intermingling backscatter behaviors. In addition, the S-1 cross-ratio (CR) shows sensitivity to phenology [43,44] and leaf loss [35]. Thus, it will be used alongside VV, VH, and Span (combined VV and VH), the dual-pol radar vegetation index (RVI) [38,45], and the Entropy-Alpha space ($H\alpha$), the latter being derived from the dual-polarimetric eigendecomposition [46].

This study explores the capability of the S-1 dual-polarimetric SAR time series specifically for assessing tree canopy damage in a deciduous broadleaf forest (DBF) under strong drought conditions in the years 2018–2020 in central Germany. The signal trajectories before, during, and in the aftermath of the drought (totaling four years) were used to investigate the physical relation of damage evolution with the SAR signal. As damage reference categories, an optically derived very high-resolution mortality dataset was created that characterizes the canopy development (CD), providing damaged and non-damaged tree sets. Finally, we display the co-evolution of the Sentinel-2 (S-2)-based NDVI and S-1 signal during the drought phases to show the complementarity of S-1 SAR to optical data for drought-induced tree canopy damage assessment.

The following hypotheses are formulated:

• A densely-sampled C-band backscatter allows for detecting tree canopy damage characterized by the mortality of individual trees in an otherwise intact canopy of a broadleaf forest.
• C-band SAR is influenced by the plant water status and structural changes in a tree canopy.
• S-1 polarimetric variables offer insight into the structural changes under such circumstances.
• Generally, C-band SAR complements measurements of optical instruments in describing temporal patterns of damage.

Two methodological implementations aim to improve the pre-processing of the dense SAR time series:

• Reducing the influence of speckle for high-resolution SAR studies via temporal-only filters.
• Implementing a novel technique for correcting the geoposition ambiguities of forest canopies in SAR geometry with a lidar surface model.

2. Materials and Methods

2.1. Study Site and Drought Impact

The study was conducted at Hainich National Park (HNP) in western Thuringia, Germany (lat: 10.4229${}^{\circ }$, lon: 51.0591${}^{\circ }$, mean elevation: 413 m a.s.l.), a mixed deciduous broadleaf forest dominated by European beech (Fagus sylvatica, 65%), ash (Fraxinus excelsior, 25%), maple (Acer pseudoplatanus, and Acer platanoides, combined 7%) with other deciduous species [47]. Due to its foundations in middle and upper Muschelkalk (shell limestones) and lower Keuper (claystone), the soils are relatively shallow, carbonated, and humus-rich with strong runoff (e.g., Rendzina, Regosol), with characteristics more pronounced on exposed slopes [48]. Thus, soils have limited storage capabilities and plant-available water rapidly declines under drought conditions.

This study refers to the successive droughts that took place in Central Europe in the years 2018–2020 [49]. Local meteorological data from an eddy covariance flux tower in the HNP indicate characteristics of an atmospheric drought and hydrological drought when compared to the 2000–2014 multi-year average [50]. Both, the anomaly of vapor pressure deficit (VPD) in 2018 of 0.302 kPa (+59%) and the extremely low SPEI drought index (description in Figure 1) indicate a strong negative water balance and atmospheric drought in these years (Figure 1). The strength of the hydrological drought was indicated by a soil water anomaly of −7.93% (−25%) in 2018, with storages not replenishing as of 2020. The climatological reference period of 2000–2014 was constrained by the measurement period of the Eddy covariance flux tower (DE-HAI), which was erected in 2000. We decided on tower data instead of a long-term record from the next meteorological station due to the large distance to the next station ( 10 km, Eisenach) and to account for DBF-specific energy fluxes at the site. We considered a 4-year period that started in 2017 and included the pre-event and the onset, as defined by the start of the SPEI anomaly, as of 4 April 2018, and the late stage of the drought (until 2020). The relative nomenclature used in this study ranged accordingly from $Y_{-1}$ (2017) to $Y_2$ (2020). It is clear that data from such a narrow study period cannot capture long-term trends and is especially prone to uncertainty. Since the Sentinel-1A/B data record is limited, we restricted the study period to 4 years.

The multi-year drought has led to substantial tree canopy damage and dieback in vast areas in Central Europe, including Germany [12,53], and was also observed in the HNP [54] (Figure 2). We confined the study site to an area of plateaus in the western HNP (Figure 3) where increased crown damage could already be successfully delineated using drone and Sentinel-2 data [54,55].

2.2. Sentinel-1 Data Processing

Only SAR data (acquired in the early morning/predawn with a local overflight time of ~7:30 CEST (6:30 CET, 5:30 UTC)) were used. At predawn conditions, healthy plants establish a hydraulic equilibrium where the water potential between the soil and leaf tissue converges (tissues replenish) ([58], p. 556). Conversely, under the effect of droughts, leaf tissues cannot fully refill during the night [59]. In this sense, water-related measurements at predawn are expected to be the most reasonable proxy for hydraulic damage in the diurnal cycle, independent of photosynthetic activity in the daytime. To avoid geometric variability among orbits due to changes in incidence angles and look directions, only scenes with average sensor incidence angles $\theta$ of 36${}^{\circ }$ over the study area were considered, corresponding to typical mid-range conditions (relative orbit 66). In total, 215 ESA Copernicus S-1 scenes with an average of 6-day revisit cycles (the same look direction, time of day, and sensor incidence angle) were processed for the period 2017–2020 [60] using single look complex (SLC) data (Figure 4). The SLC format is a precondition for performing a polarimetric eigenanalysis on the complex SAR signal and for temporal speckle filtering in radar (slant) geometry. The subsequent radar-specific processing steps were performed using the S-1 SNAP processor [61], accessed via the S-1 handler package ‘pyroSAR’ [62].

After applying the precise orbit file, all scenes were co-registered and a rolling arithmetical mean temporal speckle filter was applied pixel-wise to (1) the radiometric VV-VH intensities and (2) the complex 2-channel covariance matrix $\left[C_2\right]$ [63]. We leverage the frequent revisit time of the S-1 sensor by employing a temporal-only speckle filter to increase signal stability while keeping the maximum spatial resolution of 10 m × 10 m. Subsequently, $\left[C_2\right]$ was used for deriving $H\alpha$. We iterated a filter set of the effective number of looks (ENL) $\in \{3,7,13,19\}$ for assessing appropriate temporal filter sizes for the analysis. Despite the development of advanced time series filter methods, e.g., [64,65], filtering is usually not applicable for complex numbered data. A rolling mean operand was, therefore, used here for simplicity. Using default parameters, the radiometric data were normalized to radar brightness ${\beta }^0$ and terrain-flattened ${\gamma }_T^0$[66]. With the polarimetric data, the eigenproblem for $\left[C_2\right]$ was solved as described by Cloude [67] using custom implementations with python packages ‘numpy’ [68] and ‘numba’ [69]. Finally, all data were geocoded to a $10\times 10$ m grid (UTM 32N, WGS 84) using a lidar-derived DTM and DSM, where the latter was low-pass filtered from 1 m to 5 m to avoid ambiguities due to rough canopy surfaces and crown gaps in forests [70]. With DTM and DSM representing the minimum (no forest, DTM) and maximum (canopy height, DSM) of possible geocoding reference planes, the range of DEM influences will be judged with their extrema. By using the lidar-derived DSM instead of the default DEM settings, the geolocations of canopies shifted from the DTM ground level to the prospected canopy level. We propose that this method enables an accurate overlay of the oblique acquisition geometry of SAR with nadir-looking instruments (S-2, UAV, lidar) for forest applications. This is necessary to locate single canopies onto the SAR imagery.

To assess the structural changes to the damaged trees, the additional polarimetric indices, i.e., cross-ratio (CR), radar vegetation index (RVI) [45], and Span were derived as follows:

$$\begin{array}{cc}\hfill CR& ={\gamma }_{VV}^0/{\gamma }_{VH}^0\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill RVI& =\frac{4\ast {\gamma }_{VH}^0}{{\gamma }_{VV}^0+{\gamma }_{VH}^0}\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\gamma }_{Span}^0& ={\gamma }_{VV}^0+{\gamma }_{VH}^0\hfill \end{array}$$

(3)

where ${\gamma }_{XX}^0$ represents either VV or VH calibrated and radiometric terrain corrected S-1 backscatter.

2.3. Optical Reference Data Processing

Optical remote sensing data are sensitive to DBF crown damage in the HNP [54] and are, therefore, used as indicators for canopy damage. First, very high resolution PlanetScope imagery is used to select damaged tree canopies bi-temporally before and after the drought to create an a-posteriori canopy damage sample set. Second, the multi-spectral ESA Copernicus Sentinel-2 (S-2) time series is used alongside S-1 to understand similarities and differences between SAR and optical forest drought signals.

The PlanetScope/Dove 4-band color-infrared data were retrieved at Level-3B as an analytic ortho-scene image product at 3 m of resolution [56]. Four cloud-free scenes in peak growing seasons were chosen to create NDVI-based vitality layers for each observation year (Figure 5). Pre-event $NDV{I}_{Y_{-1}}$ was defined as undamaged. NDVI change layers were created and defined as the differences of the respective year to pre-event conditions by $\Delta NDV{I}_{Yn-Y_{-1}}=NDV{I}_{Yn}-NDV{I}_{Y_{-1}}$. As radiometric inconsistencies between PlanetScope sensors reduced the quality of time series data [71], this study restricts its use to 1-sigma quantile classification in bi-temporal (pre- and post-droughts, $Y_{-1}$ and $Y_2$) differences and visualizes damage patterns spatially. The a-posteriori Canopy Development (CD) tree sample sets were defined by:

CD −	:	Negative canopy development	$\Delta NDVI<\mu -\sigma$,
		(damaged, high mortality risk)
CD 0	:	Indifferent canopy development	$\Delta NDVI<\mu +\sigma and>\mu -\sigma$,
		(undamaged or recovered, low mortality risk)
CD +	:	Positive canopy development	$\Delta NDVI>\mu +\sigma$.
		(re-greening, disregarded here)

where $\mu$ is the mean and $\sigma$ is the standard deviation of the sample.

In the next step, the extracted CD sets were aligned to the S-1 data grid. After (bilinear) interpolation, these pixels were selected; they were members of either CD class (with 95% confidence interval). Due to the heterogeneity of the soils, forest age, height, stand density, and topography in the study area, we further employed a stratified sampling strategy to ensure that the samples were representative and classes had equal sizes. We filtered pixels with:

• Canopy density > 80%;
• Canopy height > 18 m;
• Slope angle < 10${}^{\circ }$.

The thresholds were derived from assessing the environmental variable distributions and by local expert knowledge to minimize radar-relevant topographical effects (e.g., slopes) and to pursue relative forest homogeneity, especially by excluding non-mature forest patches (canopy height and density). As horizontal density ancillary data, the “Tree Cover Density” Copernicus high-resolution layer (2018) was used [72]. Slope angle and a canopy height model were derived from lidar DTM and DSM [70]. Finally, valley features were masked out. With a total of n = 146 samples per class, 292 $10\times 10$ pixels (2.92 ha) were selected (CD + is disregarded as re-growth patterns are not the primary focus of the study).

S-2 data were retrieved via the Microsoft Planetary Computer (MPC) cloud platform [73], which was processed to Level-2A bottom-of-the-atmosphere reflectances and terrain corrected with Planet DEM 30 by the Sen2Cor processor and MPC. The normalized difference vegetation index (NDVI) was calculated at 10 m native S-2 resolution using the bands 4 (red) and 8 (nir), after [74]. The normalized difference water index (NDWI) was favored for water assessments [75]; however, the NDVI was used here due to its damage sensitivity and the higher spatial resolution of S-2 at red/NIR bands (10 m) in comparison to the SWIR bands (20 m), which should match the chosen S-1 resolution. As clouds often only obscure parts of the imagery, no global cloud filter was applied and a total of 270 scenes were downloaded. Instead, we cut NDVI data at sampling sites where NDVI < 0.25, as those values were assumed to represent clouded pixels. Then, local maxima in the time series were detected using findpeaks from the ‘gsignal’ package [76] and R statistical software [77]. The 34 scenes left (Figure 4) were then linearly interpolated to all S-2 and S-1 acquisition dates. This way, more scenes in the S-2 time series could be kept, with stable intra-seasonal signals retained. Season thresholds in this study were tailored to the local forest and species composition by using the NDVI inflection points in the spring and fall based on pre-event data. Accordingly, the growing season and full foliage cover were assumed for DOY 130 (10 May)–280 (7 October).

2.4. Time Series Analysis and Statistical Methods

At first, time series data of SAR were used to evaluate intra-seasonal and inter-annual dynamics of the damaged CD− canopies. For this, we used class-wise averaged means of the damaged tree class for S-1 indices and S-2 NDVI. In this context, different temporal speckle filters were used and their suitability is discussed.

Secondly, the state and progress of damage in tree canopies were analyzed by calculating anomalies (differences) between CD− and CD 0 tree samples. The benefits of using spatial anomalies instead of temporal ones include (1) the leveling-out of short-term influences on the radar signal (e.g., rain interception) and (2) only damage-related features are visible. For a temporal anomaly calculus, we considered a reference period based on the data since the sensor launch (from 2015) was too low.

With speckle filtering, temporal autocorrelation was introduced in the S-1 time series. Thus, to avoid bloated statistics by auto-correlated sampling, only a single S-1 scene with the largest speckle filter size $ENL=19$ was used for anomaly calculus and testing. For this, the most central scene (median of the acquisitions) in the peak leaf-on was used. Conversely, the NDVI time series did not carry out this autocorrelation. Hence, all NDVI peaks were used, as defined in Section 2.3.

Damage anomaly statistics were derived from differences in sample distributions of CD− and CD 0 as

$$\begin{array}{cc}\hfill \Delta {\mu }_{x-y}& ={\mu }_x-{\mu }_y\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \Delta {\sigma }_{x-y}& =\sqrt{{\sigma }_x^2+{\sigma }_y^2}\hfill \end{array}$$

(5)

where $\Delta {\mu }_{x-y}$ and $\Delta {\sigma }_{x-y}$ are the mean and standard deviations (SD) of differences (anomalies) between population classes x (CD−) and y (CD 0). Due to unequal variances, Welch’s t-test [77] was calculated, hypothesizing that the sample means between CD− and CD 0, were equal. Note that normal distribution and random sampling of the population are assumed.

The extracted S-1 dual-pol indices were assessed toward their information content based on the correlation analysis. Due to the variety of polarimetric indices, we used a simple Spearman rank correlation and principal component analysis (PCA) method to judge if features add information to the analysis or are redundant. Here, a PCA is a useful tool to estimate feature distance. As many polarimetric indices have been designed for full-pol configurations, the usefulness of each index needs to be assessed for the use case prior to including them in a dual-pol analysis.

3. Results

3.1. Speckle Filtering and Geocoding of Dual-Pol SAR Time Series Data

Time series signal disturbance arises from, e.g., variations in wind speed (triggering tree swaying, rain interception, or freeze/thaw events). As these short-term phenomena affect tree status mapping, their propagation to the filtered time series is of concern. In this study, we found that, by increasing the filter size, the radiometric stability for ${\gamma }_{Span}^0$ increased (Figure 6a,c) and SD decreased from ~2 dB to ~1.1 dB, for $ENL=3\mathrm{and}19$, respectively. The interpretation of plant water status changes at damaged canopies was, therefore, only valid when stable environmental and phenological conditions could be assumed. For the sensitivity analysis, we restricted the data to leaf-on-only conditions (as represented by $Y_{-1}$ NDVI seasonal inflection points). The remaining signal variances can, therefore, largely be explainable by changes in the plant water status.

The dual-pol entropy approaches a stable boundary with means of 0.7–0.8 and SD of ~0.08 at above $ENL=13$ (Figure 6b). While Lee et al. [78] suggested an $ENL=49$ for a reasonably unbiased spatial filter size for entropy, our data indicate a higher efficiency of speckle reduction in the temporal rather than in the spatial space. While the effect of combined spatiotemporal filtering has not been investigated, it holds large potential for harmonizing trade-offs. The annual trajectory of entropy suggests that scattering mechanisms simplify with the growing season and diversify again toward the fall.

Owning to the trade-off between high spatial (10 m pixels) and temporal resolution, we decided on $ENL=13$ ( 80 days for our S-1 setup) for time series depictions. Yet, for the sensitivity analysis, seasonal averages are considered, represented by an $ENL=19$ of the median scene of the summer season (see Section 2.4).

Concerning the correct detection of scattered mortality patterns, both lidar elevation models, DSM and DTM, have been used to geocode the SAR imagery. We found that the $\Delta {\gamma }_{Span}^0\mathrm{DSM}$ backscatter dropped by ~0.5 dB in the late-stage drought ($Y_2$) while $\Delta {\gamma }_{Span}^0\mathrm{DTM}$ only decreased by ~0.2 dB. As CD class differentiation is only possible with DSM-geocoded data, we infer that tracking canopy changes (by an oblique SAR geometry at a very high spatial scale) requires an elevation model that reflects the top forest layer under study rather than ground topography. Further analyses are currently being performed and documented in [79].

3.2. Temporal Signal of Drought-Stressed Broadleaf Forest

Describing the signal trajectory along the drought progression is vital to understanding anomalies induced by drought. Therefore, the time series of S-1 backscatter, CR, and S-2 NDVI at damaged tree canopies (CD−) is given in Figure 7. ${\gamma }_{VV}^0$ and ${\gamma }_{VH}^0$ have been found to decrease in the spring and rise again in the fall, complying with pre-drought conditions ($Y_{-1}$). The intra-annual changes in the plant water status can partly be disentangled in the SAR time series when carefully selecting phases of constant phenology. For example, while in $Y_0$, a near-normal leaf development can be assumed (due to wet conditions in $Y_{-1}$), the emerging drought conditions reduced vitality and chlorophyll activity from peak growing season onward (Figure 7d). Similarly, backscatter ${\gamma }_{Span}^0$ in $Y_0$ decreased toward the fall (Figure 6a and Figure 7a). When assuming foliage to be fully developed and the tree structure to be constant, the observed early summer ${\gamma }_{Span}^0$ decline can be attributed to water loss. The further decreasing ${\gamma }_{Span}^0$ in the late summer, finally, mixes with structural decay in the canopies due to, e.g., premature leaf fall and wilting. Leaf-off ${\gamma }_{VV}^0$ always remained below pre-drought conditions, indicating fundamental hydrological changes to the trees, as well as in the winter. S-2 NDVI has been found to decline in $Y_0$ toward the late summer, complying with European-wide observation of early wilting patterns [53,80]. $Y_1$ and $Y_2$ display strong vitality declines compared to the reference year.

The CR describes a magnitude with phenological phases, plateauing in summer, and increasing and decreasing in the spring and fall, respectively. In this sense, the slow slope of CR in the early summer $Y_1$ (DOY 150–200) can be related to delayed leaf development, only reaching its maximum (foliage) toward DOY ~250 (early September). NDVI-based vitality peaks earlier (DOY ~170) than CR in $Y_1$, indicating that greenness and structural development do not necessarily co-evolve in drought-affected springs but show distinct developments. As signal delays or losses cannot be observed in $Y_2$, CR eventually lacks clear interpretability in drought contexts.

3.3. The Dual-Pol SAR Information Space

The dual-polarimetric capabilities of S-1 were statistically investigated using Spearman rank correlations and a PCA. All available data of classes CD− and CD 0 were used throughout all years. ${\gamma }_{VV}^0$ and ${\gamma }_{VH}^0$ were found to have a strong positive correlation of r = 0.72 (Figure 8a). In contrast, the CR was much less related to either SAR channel with r = 0.36 (VV) and r = −0.33 (VH). The correlations between CR and both eigen-based parameters were high, at −0.50 (entropy) and −0.43 (alpha), indicating a certain alignment of information content. As the dual-pol RVI essentially contains the same information content as CR, it will be disregarded in further analysis.

Results from the PCA are depicted in Figure 8. The cumulative variance reveals that the first three axes contain the main information directions in the data (Figure 8b). When looking at the relative contribution of the variables within axis 1, the dominant information contents of ${\gamma }_{VV}^0$ and ${\gamma }_{VH}^0$ (and ${\gamma }_{Span}^0$) are visible (Figure 8c). The CR is orthogonally directed (axis 2) and similarly strong. Alpha is negatively correlated with CR, yet with a weaker signal, and the effect of entropy is low. We argue that the backscatter and CR carry most of the information content due to the orthogonality of CR to the backscatter, on the one hand, and the strong negative relation of alpha to CR and the marginal effects of entropy, on the other hand. Consequently, while ${\gamma }_{Span}^0$ is used to represent S-1 radiometry, CR represents the polarimetric information content. The eigenvalues/contributions and the eigenvectors/angles do not directly imply strength in classifying damaged forests, they only hold information for deciding the most meaningful variables.

3.4. Sensitivity and Co-Evolution of SAR and Optical Data to Drought Impact

To test the sensitivity of S-1 ${\gamma }_{Span}^0$ and CR regarding tree canopy damage under drought stress, the year-wise difference between CD− (damaged) and CD 0 (undamaged) tree samples were calculated. The S-1 backscatter of damaged trees declined moderately in $Y_1$ while significantly dropped in $Y_2$ (Figure 9), with average anomalies of −0.43 dB ± 1.36 (VH), −0.52 dB ± 1.52 (VV), and −0.50 dB ± 1.38 (Span) (Figure 10 and

Table 1. Damage anomaly (mean ± SD) for S-1 variables and S-2 NDVI during the growing season. Two sample Welch tests of p-value significance levels are given as follows: < 0.0001 = ^****, < 0.001 = ^***, < 0.01 = ^**, < 0.05 = ^·, <= 0.1 = ^., > 0.1 = non.

Index Year	Year	$\mathbf{\Delta }{\mathit{\gamma }}_{\mathit{VV}}^0$ [dB]	$\mathbf{\Delta }{\mathit{\gamma }}_{\mathit{VH}}^0$ [dB]	$\mathbf{\Delta }{\mathit{\gamma }}_{\mathit{Span}}^0$ [dB]	$\mathbf{\Delta }\mathbf{CR}$ (${\mathit{\gamma }}_{\mathit{VV}}^0$/${\mathit{\gamma }}_{\mathit{VH}}^0$) [-]	$\mathbf{\Delta }\mathbf{NDVI}$ [-]
$Y_{-1}$	2017	0.049 ± 1.432	0.094 ± 1.355	0.059 ± 1.290	−0.041 ± 1.371	0.020 ± 0.027 ****
$Y_0$	2018	0.039 ± 1.510	0.021 ± 1.365	0.038 ± 1.341	0.028 ± 1.573	−0.005 ± 0.014 ****
$Y_1$	2019	−0.182 ± 1.385	−0.135 ± 1.332	−0.172 ± 1.251 ·	−0.049 ± 1.292	−0.089 ± 0.041 ****
$Y_2$	2020	−0.524 ± 1.519 ****	−0.427 ± 1.357 ***	−0.504 ± 1.378 ****	−0.085 ± 1.563	−0.068 ± 0.031 ****

Bold entries are highlighted in the text.

). Notably, $\Delta \mathrm{NDVI}$ anomaly was largest in $Y_1$ (−0.089 ± 0.04) and decreased again in $Y_2$ (−0.068 ± 0.03), indicating that the trajectory between SAR and optical indices in late-stage drought damage diverged. It is clear from the data that the anomaly of NDVI is more certain than S-1-derived indices. However, as backscatter trends in $Y_2$ become significant, we can assert systematic changes to the scattering intensity of the C-band signal. The CD classes were defined on the basis of NDVI (PlanetScope-based) and, thus, changes to NDVI are self-evident.

The S-1 indices were compared regarding their discriminative strengths to drought damage. While $\Delta {\gamma }_{VV}^0$ had the largest anomaly and $\Delta {\gamma }_{VH}^0$ had the lowest SD among the backscatter data, $\Delta {\gamma }_{Span}^0$ had the highest t-value (4.42) in contrast to VV (4.17) and VH (3.80), indicating a higher ability for CD class differentiation. Although deviations between the signals are arguably small, we can still infer that ${\gamma }_{Span}^0$ is statistically more stable than VV and VH for detecting drought effects in this context. Regarding $\Delta \mathrm{CR}$, anomalies remained statistically insignificant over the entire observation period. Furthermore, the statistical reliability of CR anomaly was hampered here due to very high errors with SD > 1.5 ($Y_2$). Besides the potential for qualitatively interpreting changes to the seasonal cycle, CR does not show sensitivity to DBF forest droughts.

The differences between S-2 NDVI and S-1 Span signal trajectories in the drought aftermath allow for tracking multi-year damage recovery. Understanding the recovery patterns is useful in understanding the long-term implication of drought disturbances. On the NDVI/Span-plane (Figure 11), the trajectories of plant vitality versus microwave radiometry are given. Zero anomalies are expected for two cases: (1) No drought apparent, (2) drought influence for both populations (CD− and CD 0). While negative anomalies measure losses in vitality/chlorophyll activity or structural/water status for $\Delta \mathrm{NDVI}$ and $\Delta {\gamma }_{Span}^0$, respectively, positive anomalies can be caused by sampling errors or environmental unknowns. In contrast, $Y_0$ exhibits a strong drought effect for both sample populations (case 2). Only in late summer did $\Delta \mathrm{NDVI}$ decrease, which may indicate drought damage accelerating in the following years. The effects of $\Delta {\gamma }_{Span}^0$ in $Y_0$ cannot be readily interpreted because of its large fluctuation (also visible in Figure 10). As all trees are damaged (see Figure 5), the signal could be explained by random variability in tree crown desiccation as seen by C-band SAR. While NDVI shows the largest anomaly ( −0.1) in $Y_1$, ${\gamma }_{Span}^0$ anomalies are moderate ( −0.18 dB). Finally, with $Y_2$, two converse effects become visible: While anomalies in ${\gamma }_{Span}^0$ further intensify, NDVI values regain, forming a hysteresis-like circular pattern.

4. Discussion

The sensitivity of Sentinel-1 SAR for capturing tree canopy damage and mortality patterns in Central European deciduous broadleaf forests was tested. We found that the S-1 Span was the best predictor of damage detection, and the largest backscatter anomaly was found in $Y_2$ of −0.5 dB ± 1.38, indicating coupled hydrostructural decay in drought-affected trees. Polarimetric indices, such as the widely used cross-ratio (CR), could not prove sensitivity to damage within our analysis. In light of the data fusion in forest monitoring, we found the SAR anomaly time series lagging behind NDVI in time by ca. one year after the onset of the drought, suggesting complementary information content of SAR to optical sensors.

4.1. The Role of SAR Processing in Drought Observations

The applied temporal-only speckle filter scheme reduces short-term disturbances in the S-1 signal time series. A filter size of $ENL=13$ (~80 days) was shown to be a good trade-off between temporal decorrelation with phenology and biomass, and sufficient speckle reduction (Figure 6) for backscatter and polarimetric data. This enabled using credible S-1 maximum 10 m spatial pixel spacing, promoting as little spatial overlap as possible between damaged and non-damaged tree samples with their surroundings, acknowledging the heterogeneity of mortality patterns. However, with the long temporal correlation length introduced by filtering, weekly-to-monthly plant water status changes cannot be assessed independently of neighboring image acquisitions. Therefore, the signal can only be related to seasonal trends of desiccation and rehydration in canopy tissues. To increase the temporal resolution, combined spatiotemporal filtering strategies are, therefore, advised, especially in cases of high spatial homogeneity (e.g., large-scale bark beetle outbreaks). An 80-day averaging comes with the expense of ignoring short-term dynamics where outliers could corrupt the time series averaging, especially wind swaying [22], rain interception [81,82], and dew ([39], p. 508) [83,84]. To account for these factors, either threshold-based filters are found to be easy-to-use assets to gain more signal stability [82] or analytically based approaches (e.g., incorporating rain interception models [81]). Additionally, the full temporal resolution of S-1 data can be used by incorporating all relative orbit acquisitions to increase the observation density to up to 1.5 days. While this benefits the temporal coverage and allows similar ENL within a shorter time span, the incidence, azimuth angle, and diurnal variations (morning/afternoon overflight) will be introduced, which need to be further taken into account. For example, Weiß et al. [33] mitigated incidence angle influences by incorporating $\theta$ into microwave radiative transfer models.

4.2. Detecting Hydrostructural Changes to C-Band SAR in a Damage Forest Canopy

Drought conditions reduce canopy water content and limit tree leaf development and foliage duration (e.g., due to early leaf fall). Since S-1 summer anomalies of ${\gamma }_{Span}^0$ correlate with tree canopy damage, coupled hydrostructural declines in those trees are suggested (see Figure 10). We hypothesize that the observed backscatter decline can be related to (1) losses in water content in damaged plant components, e.g., due to disturbed hydraulic uptake and (2) the decline in the number of vegetation scatterers due to losses in the tree structure (leaf and, eventually, branch loss), especially in the late stage of mortality. Furthermore, since Walthert et al. [59] showed that xylem damage occurs before leaf loss in damaged mature beeches, changes in growing season backscatter may be related to mortality-induced plant water status changes instead of changes to leaf biomass, i.e., a decline in plant water content and uptake. To validate the interpretation of SAR signal losses, the link to in situ tree physiology is essential. When relating the thresholds of individual water stress to SAR, we gain the unique ability to understand plant water status and, finally, mortality in more spatial and temporal details.

As we have seen, the backscatter decline in $Y_0$ can be interpreted as the loss of canopy water content when assuming the biomass is constant (see Figure 7). The assumption is, however, violated with (early) leaf loss as the phenological boundary. With leaf fall, signal attenuation from an undamaged canopy diminishes, and contributions of the branch, stem, and soil moisture simultaneously increase [41]. For damaged trees, the winter backscatter is found to be strongly reduced (Figure 7a). The winter/leaf-off backscatter in an undamaged forest usually exceeds values measured in summer/leaf-on months. Our data suggest contrasts with normal behaviors found in the pre-drought period $Y_{-1}$ and literature (see [85]). The tree structure decline in the backscatter relates to the low water content in the desiccated damaged tissues under leaf-off conditions. Its correlation with tree desiccation is, however, conceivable. As the soil moisture of the upper layers is not variable in the winter—which may not be accurate for deeper hydrological layers—changes in soil moisture cannot explain the low winter backscatter alone. To test this, detailed knowledge about the leaf-off water content of damaged and undamaged trees, the allocation of scattered contributions to specific plant layers, and quantifying the soil contribution and soil water content are necessary. Further, physically modeling the backscatter as a function of phenology can broaden our understanding of the DBF backscatter in the winter; a crucial asset for interpreting S-1 time series data in this biome.

A promising approach exploits the separation of plant water status, phenology, and biomass dynamics when growth is absent, e.g., at a diurnal level. Here, diurnal measurements yielded encouraging results for quantifying water stress [13,15,86,87]. The suitability of S-1 for entering the diurnal observation domain is, however, far-fetched due to the long time lags between morning and evening acquisitions (3–6 days). Approaches for forests, while successfully used for water stress detection in agricultural areas [88], need to be developed. As mentioned earlier, using different S-1 orbits further introduces variabilities in incidence and azimuth angles that need to be considered.

4.3. Potential of S-1 Polarimetry for Detecting Changes in Scattering Mechanisms

Changes in scattering mechanisms in drought-damaged tree canopies were expected to shift the polarimetric SAR signatures from volume scattering to simpler scattering mechanisms after the onset of drought. Dual-polarimetric S-1 variables (CR, RVI, entropy, alpha) showed sensitivity to the seasonal dynamics of deciduous trees (see entropy in Figure 6b and CR in Figure 7c) and, thus, covaries with phenology. This supports the interpretation that CR is sensitive to the presence and absence of leaves and, thus, reacts to substantial changes to the tree structure [32,43,85]. Still, changes in CR with crown damage could not be observed in our study. In contrast to the findings in [35,36], which used CR as a proxy for defoliation, similar signs of signal loss were not significant in our data during the growing season (see CR anomaly in Table 1). We conclude that CR helps to interpret phenology in DBF forests. However, a clear role in the structural decay of CR under drought stress needs further research and detailed knowledge about the canopy structure and phenology.

Due to the high-entropy regime (H > 0.6), ambiguities in scattering mechanisms are common problems for forests, especially at the C-band, impeding the successful retrieval of polarimetric information. Our experimental results comply with the simulations by Ji and Wu [89], showing that co-cross (VV-VH) polarization combinations could not correctly identify the canonical scattering mechanisms compared to full-polarimetric or co-pol dual-polarimetric setups (VV-HH). Still, a careful 3D separation of backscatter contributions from tree compartments (twigs, branches, stems) and the ground would strongly facilitate interpreting the detected annual cycle of the dual-pol entropy. Studying forests with multi-temporal, fully polarimetric SAR systems would help researchers understand microwave scattering in forests under changing phenology.

CR and Span were statistically found to be the most independent in the PCA based on the orthogonality of information content. As either variable can be directly calculated from VV and VH intensities (e.g., from ESAs GRD product), one can reconsider sparring the computationally demanding processing of complex-numbered SLCs and the eigenanalysis in favor of the simpler CR and Span. The high magnitudes and low errors found in Span indicate a slight improvement over using VV or VH channels alone in terms of more stable anomaly detection. Therefore, we suggest its use for mortality assessment and, more generally, for applications focusing on SAR radiometry.

4.4. Time-Lagged Damage Patterns in the NDVI/Span-Plane

An NDVI/Span-plane was introduced that reveals the lagged signal anomaly in the backscatter against optically-derived vitality indices along the impact and recovery phases of the drought. If trees desiccate and die, water content and leaf/twig loss will lead to long-term structural decay, as captured by the S-1 Span. At the same time, photosynthetic activity showed tendencies for recovery ($Y_2$). Thus, a time-lagged correlation between $\Delta \mathrm{NDVI}$ and $\Delta {\gamma }_{Span}^0$ is implied, highlighting the different aspects and capabilities of optical and SAR sensors for detecting damage patterns and, ultimately mortality. While this approach is not suitable for characterizing drought sequences, it can be the basis for the multi-sensor remote sensing of disturbance ecology.

The NDVI/Span-plane can help researchers understand coupled photosynthetic and hydrostructural decay in the time domain for studying the damage and recovery of forests. In this sense, a forest would return to its Cartesian origin once it recovers from drought disturbance. The position of the hysteresis would then characterize its drought impact and recovery stage.

The recovery of $\Delta \mathrm{NDVI}$ in $Y_2$ suggests an increase of greenness and vegetation vitality at damaged tree sites (CD−) at the late stages of droughts. A possible explanation is the development of undergrowth vegetation and pioneering plants underneath desiccated/dead trees. Moreover, if droughts are not lethal for individual trees, the partial regrowth of compartments may occur, or mixed pixels containing adjacent recovered or undamaged individuals may affect NDVI. Positive values of anomalies between CD− and CD 0 are found in $Y_{-1}$ and $Y_0$, indicating inconsistencies in the sample set. The data cannot explain positive values when assuming that the samples are drawn from the same tree population and damage leads to a decline in the backscatter and NDVI. Unknown variations in the population can play roles, including demography (age, species), forest structure (vertical and horizontal distribution), topography (e.g., aspect), soil/geology (e.g., soil composition, rock formations, runoff scheme), and variables that go beyond the initial sampling scheme (NDVI threshold, slope, canopy height, tree cover density).

5. Conclusions

Our results suggest that the Sentinel-1 backscatter can track the evolution of tree canopy damage at very high spatial (10 m × 10 m) and moderate temporal (~80 days) scales in a temperate deciduous broadleaf forest. S-1 polarimetric indices did not show the significant separation of damaged and undamaged trees, which was partly attributed to the high entropy scattering in the forests and the dual-pol sensor properties. The S-1 Span showed the best capability for assessing tree canopy damage. If considering periods of constant biomass, such as the peak growing season, backscatter anomalies can be carefully related to the plant water status. However, additional data on the plant water status are necessary to validate this conclusion. With the SAR anomaly time series lagging behind optically derived vitality indices in time by ca. one year, the long-term structural decay in the aftermath of drought years can be monitored. Further, with the help of the circular pattern in a Span/NDVI observation space, a tool for quantifying forest decline and recovery stages is provided. Based on the introduced methodology, a SAR-based or SAR-enhanced monitoring tool for observing long-term forest drought damage and its application in other broadleaf is conceivable. For mapping damage, integrating the sensors would enable the unique benefits of each frequency domain for retrieving mortality-relevant forest parameters, by, e.g., combining NDVI-based vitality and phenology information with SAR-based plant water status and structure.

Since speckle variance strongly decreases confidence in anomaly detection and still hinders out-of-the-box mapping with our data, new methods for estimating speckle patterns in time should be applied when working at very high spatial resolutions. Extending the reference dataset to other forest types and areas of known tree mortality would build the basis for SAR damage assessment outside of this local experiment. Other central unknowns comprise the discrimination of plant water status from meteorological factors, tree structure, and phenology, and understanding the backscatter contributions from different vegetation layers in the C-band, primarily under drought stress. At this point, the data gap in a tree’s physiological measurements for retrieved remote sensing forest hydrological parameters strongly limits advancing toward forest stress mapping. Shortcomings in the interpretability of SAR signals should further motivate research to analytically model microwave scattering along phenological and diurnal phases to retrieve the stand-level water status directly from microwave measurements. While recent advances in VWC modeling based on the transmission parameter VOD show promising results for linking water status and drought stress to microwave measurements at a coarse spatial resolution, the transfer to active SAR systems will be a crucial advantage for high-resolution mapping of forest health and tree mortality.