Satellite data reveal differential responses of Swiss forests to unprecedented 2018 drought

Abstract Extreme events such as the summer drought of 2018 in Central Europe are projected to occur more frequently in the future and may cause major damages including increased tree mortality and negative impacts on forest ecosystem services. Here, we quantify the response of >1 million forest pixels of 10 × 10 m across Switzerland to the 2018 drought in terms of resistance, recovery, and resilience. We used the Normalized Difference Water Index (NDWI) derived from Sentinel‐2 satellite data as a proxy for canopy water content and analyzed its relative change. We calculated NDWI change between the 2017 pre‐drought and 2018 drought years (indicating resistance), 2018 and the 2019 post‐drought (indicating recovery), and between 2017–2019 (indicating resilience). Analyzing the data from this large natural experiment, we found that for 4.3% of the Swiss forest the NDWI declined between 2017 and 2018, indicating areas with low resistance of the forest canopy to drought effects. While roughly 50% of this area recovered, in 2.7% of the forested area NDWI continued to decline from 2018 to 2019, suggesting prolonged negative effects or delayed damage. We found differential forest responses to drought associated with site topographic characteristics and forest stand characteristics, and to a lesser extent with climatic conditions and interactions between these drivers. Low drought resistance and high recovery were most prominent at forest edges, but also on south‐facing slopes and lower elevations. Tree functional type was the most important driver of drought resilience, with most of the damage in stands with high conifer abundance. Our results demonstrate the suitability of satellite‐based quantification of drought‐induced forest damage at high spatial resolution across large areas. Such information is important to predict how local site characteristics may impact forest vulnerability to future extreme events and help in the search for appropriate adaptation strategies.

Impacts of extreme events significantly affect ecosystems and humans at large (IPCC, 2021) through, for example, increased wildfires and air pollution, agricultural loss or green and blue water scarcity (Berdanier & Clark, 2018;Clark et al., 2016;Miralles et al., 2019;Puletti et al., 2019). Generally, the vulnerability of ecosystems increases when facing droughts (Allen et al., 2010(Allen et al., , 2015IPCC, 2014;Millar & Stephenson, 2015). However, it is not yet fully clear what effects these extreme events may have on ecosystems such as temperate forests in Switzerland and Central Europe. The extreme droughts of the 21st century (i.e., 2003, 2015, and 2018) are within the variability of drought events (Ionita et al., 2021) regarding precipitation deficits but unprecedented warming in the last millennium has led to an increase in frequency and magnitude of the compound climate extreme "hotter droughts" with disproportionate effects on ecosystems (Zscheischler & Seneviratne, 2017). After drought events, reduced primary productivity (Ciais et al., 2005), higher vulnerability to tree-mortality, forest decline (Allen et al., 2010(Allen et al., , 2015, shifts in species composition, and subsequent detrimental effects on ecosystem services  have been observed to be dependent on site conditions (Lévesque et al., 2013;Rita et al., 2020). Nevertheless, it is still unclear what the long-term consequences in terms of productivity and the tree's ability to resist pests and pathogens are (Ciais et al., 2005) and how these consequences influence global forest health (Allen et al., 2015). Furthermore, the need to extend from local field measurements to larger spatial scales (Lévesque et al., 2013) while maintaining fine-scale spatial heterogeneity (Rita et al., 2020) for more generalizable results has been asserted.
Droughts, especially in combination with increased temperatures, represent a typical case of extreme events that become more intense stressors for temperate ecosystems and are likely to occur more frequently in the future (Horton et al., 2015;IPCC, 2014;Kornhuber et al., 2019;MeteoSchweiz, 2018). Drought types, that is, meteorological (Wilhite & Glantz, 1985), hydrological (Wilhite, 2000), agricultural (Dai, 2011), and ecological (Crausbay et al., 2017), may depend on different processes and cause different responses, making it necessary to clearly define which categorization is being used for a given drought event (Slette et al., 2019;Zang et al., 2020). Here, we use the concept of ecological drought for our quantification of its impacts. Ecological drought has been defined as "an episodic deficit in water availability that drives ecosystems beyond thresholds of vulnerability, impacts ecosystem services, and triggers feedbacks in natural and/or human systems" (Crausbay et al., 2017). Water availability deficits are driven partly by meteorological drought events or precipitation deficits that decrease soil moisture, runoff, and water storage. Further factors such as temperature, wind speed, radiation, and humidity alter the atmospheric evaporative demand (Miralles et al., 2019). Therefore water availability deficit affects plants directly by diminishing hydraulic conductance (Breshears et al., 2013) and indirectly through increased evapotranspiration (Teuling et al., 2013).
These drought, responses include changed species abundances, interaction frequencies, and herbivory levels, that may even result in species turnover (Jiang et al., 2019;Sippel et al., 2018), or changed ecosystem functions such as primary productivity, erosion control, and fire suppression (Reichstein et al., 2013).
Climatic and location-specific factors influencing forest vulnerability to drought include higher temperatures, which increase the atmospheric evapotranspirative demand and thus increase water stress (Allen et al., 2015;Buras et al., 2020). Forests in more exposed locations, for example, in ridges, are confronted with generally drier soils (Hawthorne & Miniat, 2018), often driven by higher wind and radiation exposure. Fragmentation plays an important role too, as forest edges may be more vulnerable to drought effects. At edges, both the micro-climatic conditions and the larger canopy volumes result in higher evapotranspiration Buras et al., 2018). The species-specific water uptake depths of broadleaf trees (e.g., European beech, Fagus sylvatica) in contrast to conifers (e.g., Norway spruce, Picea abies) can shift to deeper soil levels and thus depleting water availability (Brinkmann et al., 2019). Higher forest density increases competition for water resources (Crowther et al., 2015;Sohn et al., 2016), whereas higher biodiversity allows functional compensation to stabilize forest ecosystems and increase resistance to disturbances (Anderegg et al., 2018;Isbell et al., 2015;Loreau et al., 2003).
Here, we aim at quantifying the effect of the abovementioned environmental influence factors on the drought response of Swiss forest ecosystems. With this, we respond to the need for comprehensive assessments to contribute to understanding whether our forests are resilient to anticipated climate change. More specifically, using Switzerland as a model system that provides many environmental gradients (Eilmann et al., 2009;Jolly et al., 2005;Rigling et al., 2013;Wolf et al., 2013), (i) we provide a remote sensing based assessment of drought effects based on canopy water content and (ii) we test the effect of the interplay of abiotic and biotic environmental drivers on forest drought responses. We took advantage of the naturally occurring sequence of pre-drought, drought, and post-drought conditions between 2017 and 2019 to assess the resistance, that is, capability to withstand disturbance (2017-2018), recovery, that is, capacity to reclaim performance after the disturbance (2018-2019), and resilience, that is, capacity to regain predisturbance levels (2017-2019) of forest ecosystems to droughts (Albrich et al., 2020;Ingrisch & Bahn, 2018;Isbell et al., 2015;Lloret et al., 2011;Wagg et al., 2017). The combination of low precipitation and high temperature during the summer 2018 caused early leaf and needle discoloration, premature leaf shedding, canopy dieback and tree mortality among many forests in Switzerland (Denzler, 2018(Denzler, , 2019Medienmitteilung Basel Landschaft, 2018;Schuldt et al., 2020;SDA, 2019). The most affected tree species included Norway spruce and European beech (Schuldt et al., 2020), the two most common tree species in Switzerland occurring across 43.7% and 18.1% of the Swiss forest area, respectively. Previous studies focusing on the 2018 extreme drought event in Switzerland, compared it with the lower temperature drought in 1947 (Rathgeb et al., 2020), examined the geographical variability of its effects (Zappa et al., 2019), addressed its effects on the dominant tree species European beech (Baltensweiler et al., 2020) and forest health (Rohner et al., 2021). Across Central Europe, satellite-based measurements of the Normalized Difference Vegetation Index have been used to assess impacts of the 2018 drought event Buras et al., 2020;Schuldt et al., 2020). But we still lack knowledge on drought responses over large spatial extents while maintaining fine-scale heterogeneity. In the present study we used Sentinel-2 satellite data to derive a time series of the Normalized Difference Water Index (NDWI) as indicator of canopy water content (Gao, 1996).
We tested whether gradients of temperature, topography, fragmentation, and species heterogeneity affected forest drought responses. The underlying hypothesis is that more drought damages as measured by lower NDWI should occur at higher temperatures because of its effect on the evaporative demand, at ridges and along forest edges because of the effect of wind and radiation, for coniferous tree species (e.g., Norway spruce), and in mono-specific or low-species heterogeneity stands as more biodiversity is expected to provide more resilience. Together, we expect that tests of these hypotheses will further elucidate drought responses in temperate forests and increase our understanding of how such ecosystems respond to this increasingly more frequent stressor.

| A natural experiment-the 2018 summer drought in Switzerland
The summers of 2017, 2018, and 2019 were among the five warmest summers on record in Switzerland since measurements started in 1864 (Bastos et al., 2020;Brun et al., 2020;Buras et al., 2020;Gharun et al., 2020;Rohner et al., 2021;Schuldt et al., 2020)  To capture the effects of the 2018 drought year on vegetation, we used the climatic water balance (CWB) as a measure of drought stress. This measure combines precipitation and temperature and, therefore, relates to water availability for plants (Senf et al., 2020;Vicente-Serrano et al., 2010. We calculated CWB by subtracting the potential evapotranspiration, derived from mean monthly minimum and maximum temperatures (Hargreaves, 1994), from the monthly precipitation sum (Thornthwaite, 1984). We used monthly precipitation and temperature data products from MeteoSwiss (named RhiresM, TminM, and TmaxM in the MeteoSwiss data catalogue) at 2 km spatial resolution (Federal Office of Meteorology and Climatology, MeteoSwiss; https://www.meteo swiss.admin.ch/ home/clima te/swiss -clima te-in-detai l/raeum liche -klima analy sen. html). We used the CWB values integrated from March to August for 2017, 2018, and 2019 to calculate CWB anomalies defined as the difference from the norm period 1981-2010 in terms of standard deviations of the norm period (z-score). For example, the CWB anomalies for 2018 were more than two standard deviations below the norm period over large areas of Switzerland ( Figure 1).

| Using Sentinel-2 data to assess forest condition
Copernicus Sentinel-2 satellite data provide radiometric information of the visible to shortwave infrared (SWIR) in 13 distinct spectral bands with a spatial resolution of 10, 20, and 60 m (Drusch et al., 2012).
We used three image mosaics comprising 11 Sentinel-2 level-1C (top-of-atmosphere reflectance) tiles to assess forest state in August 2017, 2018, and 2019. August was chosen since at this time of year there is a larger contrast between symptoms of environmental stress and the absence of senescence following the natural phenological cycle of Swiss forests (Zubler et al., 2014). Furthermore, August represents the peak time of the combined heat wave and drought in 2018 (MeteoSchweiz, 2018).
We applied an atmospheric correction to retrieve bottom-ofatmosphere reflectance (level-2A) using the Sen2Cor processor (Louis et al., 2016) of the Sentinel Application Platform (SNAP -ESA Sentinel Application Platform v6.0, http://step.esa.int). We also applied an empirical line-based radiometric co-registration of images (Ariza et al., 2018) to compensate possible inaccuracies in obtained reflectance values due to imperfect atmospheric characterizations over time and, thus, to further increase consistency of the time series.
This procedure inter-calibrates image radiometric information of a specific image of the time series considering a pre-defined reference image (e.g., the first image of the time series). Therefore, a regression line is derived for corresponding reflectance values in two images (i.e., specific and reference image) that represent invariant targets, that is, materials whose reflectance is invariant or varies only slightly across time (e.g., artificial surfaces such as concrete or asphalt). The regression line is then used to adjust the reflectance values of the entire specific image considering the reference image. Therefore, we selected 11 to 29 pseudo-invariant features for each tile (total 178 invariant locations) to estimate and correct wavelength-dependent reflectance differences between years.
We masked clouds and their shadows using simple band thresholds determined by expert knowledge on pre-processed reflectance data. Because clouds largely reflect light in the blue wavelength region, they are well distinguishable from low reflective forest areas in this wavelength region. The near infrared is well suited to distinguish cloud shadows from vegetation because vegetation shows substantially higher reflectance compared to cloud shadows. We then iteratively tested increments of 1% in reflectance in the blue and near infrared over different forest types to determine thresholds for clouds and cloud shadows. We disregarded all pixels with more than 5% reflectance in the blue (band 2 of Sentinel-2, centered around 492 nm, bandwidth 66 nm, 10 m spatial resolution) and less than 15% reflectance in the near infrared (band 8 of Sentinel-2, centered around 833 nm, bandwidth 106 nm, 10 m spatial resolution) for every single pre-processed image. It should be noted that the applied cloud shadow threshold yields some false positives, for instance, vegetated pixels may be affected by cast shadow along mountain ridges. Besides clouds, artificial atmospheric effects (i.e., contrails) and other sensor errors generating false spectra can affect the analysis. Assuming that such errors only occur temporarily and manifested forest drought responses do not change substantially within one month, we calculated per pixel and spectral band the median from all pre-processed and masked images acquired in August of a given year (1-7 images) to remove these effects and to work with one filtered spectral data cube per year.
The reflectance data obtained with the above procedure were used to calculate NDWI values (Gao, 1996) in the following way: where (λ) stands for the top of canopy reflectance at a specific wavelength (λ; nm), which differ from the original definition using 860 nm and 1240 nm instead. The use of near infrared band 8 (centered around 833 nm, bandwidth 106 nm, 10 m spatial resolution) instead of band 8A (centered around 864 nm, bandwidth 21 nm, 20 m spatial resolution) comes at a cost of lower spectral resolution but has the advantages of a higher spatial resolution, and it has been successfully demonstrated that the such calculated NDWI is sensitive to water stress (Marusig et al., 2020). Higher spatial resolution also justifies the prioritization of band 11 (centered around 1614 nm, bandwidth 91 nm, 20 m spatial resolution resampled to 10 m) over band 10 (centered around 1374 nm, bandwidth 31 nm, 60 m spatial resolution). Furthermore, the whole SWIR is dominated by water absorption (Hill et al., 2019) and although the 1614 nm feature has much lower absolute reflectance values than the 1240 nm feature, it will only make the normalized difference more pronounced in comparison to the original NDWI. This adjusted NDWI was proposed by the European Drought Observatory and is sometimes called NDWI 1640 (Wang et al., 2019), normalized difference infrared index (Wilson & Norman, 2018), or normalized difference moisture index (Jin & Sader, 2005). insufficient compensation of illumination geometries and atmospheric disturbances. In particular the actual irradiance over shaded hillslopes is difficult to estimate, resulting in wavelength-dependent uncertainties that increase with shorter wavelengths, and propagate to retrieved bottom-of-atmosphere reflectance and derived vegetation indices (Damm et al., 2015). The identification of a robust index less sensitive to such topographical effects is essential to avoid misleading results. We compared the sensitivity of the NDWI to detect drought by comparing it to several indices commonly proposed to detect plant water stress, including the green-red vegetation index (Tucker, 1979), the red chromatic coordinates (Gillespie et al., 1987), the normalized difference vegetation index (Tucker, 1979), and the enhanced vegetation index (Huete et al., 2002). We calculated these indices from the Sentinel-2 imagery for August 2018 and compared mean and variance of the index values for a selection of forest pixels divided into three categories. These categories were as follows: (i) forest appearing healthy on flat terrain, (ii) forest visibly damaged (leaf discoloration) on flat terrain, and (iii) forest appearing healthy on steep slopes. A good index that would not be sensitive to topography should have similar values for the two categories (i) and (iii) but not (ii). This condition was best met by NDWI ( Figure S1), likely because NDWI involves wavelengths in the near and short-wave infrared that are less affected by wavelength dependent uncertainties common in the visible region (Guyot et al., 1992).

| NDWI change from pre-to post-drought as measure of forest drought response
To assess the response of forests to the 2018 drought year, we calculated NDWI changes (ratios) between (i) 2017 and 2018 as measure of resistance, (ii) 2018 and 2019 as measure of recovery, and (iii) between 2017 and 2019 as measure of resilience (Albrich et al., 2020;Ingrisch & Bahn, 2018;Isbell et al., 2015;Lloret et al., 2011;Wagg et al., 2017). Preceding dry conditions affecting the base condition (defined as August 2017) or recovery lasting longer than a year either from the drought 2018 or previous events were disregarded.
To avoid autocorrelation and to aggregate the very large number of pixel data for further analysis, we then converted the continuous change measures into binary data using a 10% change as cutoff between "no change" and "strong change," distinguishing between negative and positive changes. This binary approach divides pixels into non-damaged and damaged forest, a differentiation important for our application purposes. The threshold allowed us to count the number of changed pixels relative to all pixels within a stratum containing pixels from a certain area or falling into the same category of an environmental variable (e.g., distance from forest edge). The percentage of positively or negatively changed pixels within a stratum was then used as the dependent variable measuring the response of forest to drought. The cutoff of 10% was chosen in the following way. We first flagged all forest pixels with leaf discoloration visible in high-resolution orthophotos from August 2018 for a spatial subregion of ~10 km 2 (Geoinformation Kanton Zürich, Orthofoto Sommer RGB/Infrarot 2018, https://maps.zh.ch/). Then, we iteratively tested increments of 1% change in NDWI until the resulting maps would cover the full extent of damage identified in the orthophotos. no change would be indicative of this response; however, using our measure, "over"-resilience indicates that positive recovery "overcompensates" negative resistance whereas "under"-resilience indicates that forest condition was still reduced in the post-drought year compared with the pre-drought year.

| Validation
We validated the accuracy of the detected tree damage with highresolution satellite imagery acquired around the 2018 drought event and retrieved from Google Earth Pro (https://www.google.com/intl/ en/earth/). We selected sample sites at regular intervals (~30 km) although missing summer imagery around the drought event resulted in gaps. Particularly in the alpine areas, imagery was scarce, whereas the northern and western part of Switzerland were widely covered.
We identified damage in tree crowns that are related to canopy water content and thus NDWI (e.g., leaf discoloration for broadleaf trees ( Figure 2b) or needle shedding for conifer trees (Figure 2e)) and digitized the extent of affected tree crowns. We identified and digitized tree crowns with no visible damage in close vicinity (max 4 km).
In result, we collected 69 of these pairs of damaged and healthy trees for validation purposes. The accuracy of detected tree damage was assessed by the overall accuracy, the producer's accuracy, and user's accuracy from a standard confusion matrix.

| Environmental variables potentially influencing drought response of forests
Switzerland offers large gradients of environmental conditions that may impact responses of forests to drought .
We therefore assessed environmental variables that could potentially be the most important drivers of forest responses to drought.
First of all, we assessed drought stress severity itself using climatic environmental variables. We used temperature and precipitation from March to August 2018 (Federal Office of Meteorology and Climatology, MeteoSwiss; https://www.meteo swiss.admin.ch/ home/clima te/swiss -clima te-in-detai l/raeum liche -klima analy sen. html), as well as the above-defined CWB anomaly (z-score) (see section A natural experiment-the 2018 summer drought in Switzerland) to assess hotter (higher temperatures), drier (lower precipitation), and hotter and drier (low CWB z-score) conditions. We aggregated monthly data products of the growing season (March-August) for mean temperature and total precipitation of the drought year (named TabsM and RhiresM in the MeteoSwiss data catalogue, 2 km grid size, resampled to fit Sentinel-2 pixels) and compared the aggregated values with their corresponding values from 1981-2010 as percentual differences.
Next, we used data from the swissALTI 3D digital elevation model (DEM) (Federal Office of Topography, Swisstopo; www.swiss topo. admin.ch/en/geoda ta/heigh t/alti3d.html) to derive elevation, slope, and folded aspect for each forested pixel of 10 × 10 m. We used the folded aspect to quantify to which degree slopes were facing south.
The slope, aspect, and latitude were used to calculate the potential direct incident radiation (PDIR) (McCune, 2007;McCune & Keon, 2002). Additionally, we derived the topographic position index with an annulus (inner radius 150 m and outer radius 300 m) from the DEM (Weiss, 2001), where positive values indicate higher exposition.
We used the Swiss National Forest Inventory (NFI) forest mask with a spatial resolution of 1 x 1 m (Waser et al., 2015;WSL, 2020) to distinguish forested from non-forested areas. The total forested area of Switzerland is 31% (BAFU, 2019). We applied a 20 m inward buffer to the forest mask to exclude mixed pixels due to the coarser spatial resolution of the satellite data compared with the NFI data and possible geometric distortions. We also used the forest mask to determine the closest distance of a forest pixel to the forest edge.
We used the NFI tree type map (Waser & Ginzler, 2018;Waser et al., 2017) to distinguish broadleaf and coniferous trees and approximate species heterogeneity within 100 × 100 m moving neighborhoods.
Overall, 28.6% of all Sentinel-2 pixels covering Switzerland, not excluded due to clouds, topography, or buffering of the forest mask, represent forest and were further analyzed.

| Statistical analysis
First, we used simple linear mixed models with ordinary least squares to assess the relationships between percentage affected forest pixels (resistance, recovery and resilience, for positive and negative changes; values were angular transformed [arcsine of the square root of the proportion] to ensure normality and homoscedasticity of residuals) and single environmental variables (Table 1) Figures S2-S4). In addition to analyzing the percentages of pixels with >10% change between years, we also analyzed the mean change, that is, mean resistance, mean recovery, and mean resilience directly ( Figure S5), and to verify the used 10% threshold we also tested 5% ( Figure S6) and 15% ( Figure S7) cutoff values. Because all these analyses yielded similar results, we focus our presentation on the analysis with 10% cutoff values.
To obtain a conditional analysis for recovery, we repeated the above analyses for the subsets of pixels that had only positive, only negative, or only no changes >10% from 2017 to 2018. This allowed us to test if, for example, positive recovery was generally following low resistance, that is, affecting pixels that had lost >10% in NDWI from 2017 to 2018, or if frequently pixels would be negatively affected in the two yearly intervals in a row.
Second, we used more complex linear mixed-effects models to assess the effect of the most significant environmental variables from the univariate analysis described above as well as their interactions. However, to avoid overfitting with too many strata (bins), we made an aggregated multidimensional table using three intervals for each of five chosen environmental variables and the 10 regions (see description below). With this reduction of the dataset from millions of pixels to 3 5 * 10 = 2430 strata or bins-of which 2108 were non-empty and for which the percentage affected pixels could be calculated-the data could be handled by lme() to test significances and lm() to obtain effect sizes. For the latter we used percent SS, which correspond to increments of multiple R 2 as terms are fitted in sequence. Note that changing the sequence of terms had very little influence because the aggregated table resulted in high orthogonality among explanatory variables even when using number of pixels per stratum or bin as weighting variable. We used three intervals per environmental variable to allow for detection of potential non-linearities and, in cases where we detected these, we accounted for them using contrasts . The intervals for the environmental variables were the following: for elevation:   Note that when a contrast term is followed by the factor, the latter provides a test for remaining differences among levels once the differences due to the contrast have been explained. For example, ele + ELE partitions the factor ELE (2 degrees of freedom (DF)) into the contrast lowland versus foothills and Alps (ele, 1 DF) and the remaining difference foothill versus Alps (ELE fitted after ele, 1 DF). Thus, the models (i)-(iii) all fit the same number of parameters, namely 51 (see Table S1). Fitting main effects and all 2-way interactions for the five environmental terms, we could keep the number of fitted parameters small relative to the number of data values (more than 40 data points per estimated parameter, that is, DF of residual 2048). As for the single analyses above, we also compared the results of the complex analysis using percentages of pixels with >10% changes between years with a complex analysis with the mean changes in NDWI between years directly, that is, mean resistance, mean recovery, and mean resilience (Table S2).  High NDWI increase (%) (6.86%, Figure 3a). The lowest proportion (1.63%) was found in the South. For 3.76% of the forest area, NDWI values increased by ≥10%.
The highest proportion of areas with increasing NDWI from 2017-2018 (here referred to as "over"-resistance) were found in the southeastern part of Switzerland (11%, Figure 3d). Generally, 2.15% of forest areas showed strongly increasing NDWI from 2017 to 2019 (here referred to as "over"-resilience) with 3.37% along the northern side of the Alps. This suggests that the temperature component of drought can actually be beneficial in areas with generally high humidity and low temperature (e.g., north-facing slopes at high elevation).

| Relation between further environmental variables and forest drought response
Compared with the above climatic variables, topographic variables and variables related to forest structure showed stronger linear relationships with the proportion of forested area with strong NDWI changes between years ( Figure 4). The proportion of forested area with ≥10%NDWI decreases from 2017-2018, that is, low resistance, was strongly increasing from northern to southern aspects ( Figure S2) and towards forest edges (for linear model despite a decreasing decline with larger distances as seen in Figure S3). The proportion of forested area with ≥10% NDWI increases from 2017-2018, that is, over-resistance, increased along with elevation, was large on northern aspects and in forests with a large proportion of coniferous trees or higher species homogeneity. As mentioned in the previous section, these are likely areas with generally high humidity and low temperature and thus might have suffered less from the lowered precipitation and actually benefitted from the increased temperature component of the drought event.
The proportion of forested area with ≥10% NDWI decreases from 2018-2019, that is, negative recovery, was strongly increasing with the proportion of coniferous trees but also with elevation, steepness of slopes, and exposed positions. However, positive recovery was positively related to forest edge proximity, southern aspect, and high radiation, that is forest areas that were more negatively affected by the drought in the first place (low resistance, see above; in the next section recovery is analyzed separately for these areas with low resistance).
Regarding resilience, where no strong changes in NDWI from 2017 to 2019 indicate compensating resistance and recovery processes or generally small NDWI changes between adjacent years, strong changes  Table 1 (max-min rather than min-max) to focus on drought characteristics (climate) or site characteristics where higher proportions of damaged forest were found (forest stand). The intensity of color represents the effect size as measured by the R 2 value of the regression. Fields of non-significant relations (p > .05) are empty. Corresponding analyses with mean NDWI changes and with NDWI changes of ≥5% or ≥15% are presented in Figures S5-S7, respectively. These additional analyses yielded similar results to the ones presented here strong NDWI increase from 2017 to 2019 (over-resilience) was large at forest edges (R 2 = .59), in higher elevations, less exposed positions and again in forests with a high proportion of coniferous trees ( Figure S4).
The results from the ≥10%-change analyses presented in the previous three sections were further supported by the results from the direct NDWI-change analyses ( Figure S5) and from analyses using ≥5% ( Figure S6) or ≥15% ( Figure S7) changes in NDWI. near the forest edge, on south-facing slopes, and for broadleaf trees ( Figure 5). We generally found less recovering forest at higher elevation, on steeper slopes, at locations with lower potential incident radiation and in more homogenous forests, although these relationships had low R 2 values. The relationships were all inverted for forest areas that either experienced further NDWI decrease (NDWI decrease 2018-2019 ≥10%) or did not change by more than 10% in either direction from the lowered value after the drought event.

|
Thus, areas with no or even negative recovery from 2018 to 2019 were found on north-oriented slopes, at more exposed locations, further away from the forest edge, and in forests with a greater proportion of coniferous trees.
Proportions of forest areas with either negative (2%) or positive (3%) recovery were small for forests with stable NDWI values from 2017 to 2018 ( Figure S8). However, forest areas with ≥10% increased NDWI values from 2017 to 2018 often had negative recovery (23%), indicating compensation in the opposite direction from the typical one reported above of low resistance followed by high (= positive) recovery.

| Combined effects of environmental drivers and their interactions
The analyses of the combined effects of the strongest environmental drivers from the separate regressions presented above confirmed the pattern found there, although here each variable was only represented with three intervals, here defined as categories, along its entire range.
This method-which can be used for bias reduction in observational studies (Imbens & Rubin, 2015)-had the advantage that correlations between environmental drivers could be largely removed and thus the fitting sequence of drivers had little influence on their contribution to multiple R 2 values, which we used to assess the contribution of single environmental drivers or their interactions to forest responses ( Figures 6 and 7). Thus, all tested environmental drivers (i.e., elevation, folded aspect, topographic position index, distance to the forest edge, and forest type), but not so much their two-way interactions, were generally highly significant in linear mixed-effects models with the factor region as random-effects term (Table S1). The combined analyses had multiple R 2 values between .73 and .88, which is high compared with the large number of strata used in the analysis (see section Statistical analysis). A large contribution came from the random factor region (fitted first when calculating % SS), which was expected due to its higher degrees of freedom, but even after accounting for these regional differences, contributions from the tested environmental drivers were still large, together explaining typically around 40% (R 2 * 100-variance explained by regions in Figures 6 and 7).
The combined analysis showed that the forest structure variables, distance to forest edge and proportion of coniferous trees were the most important drivers of forest responses to the 2018 drought; among the topographic variables, aspect and to a lesser degree elevation were important drivers. Using contrasts, we found that for distance to forest edge the largest differences in forest response was between forest edge (≤20 m) and forest interior. For forest composition, the largest differences were found between broadleaved and other forests, and for aspect between southern and other directions. However, elevation effects were only partly explained by the contrast lowland versus foothills and Alps, sometimes leaving larger differences between foothills and Alps (right panel in Figure 7).
There were only three relatively strong interactions: forest edges in lowland areas (interaction "ele:dis") had particularly low resistance (left panel in Figure 6) but high recovery (middle panel in Figure 7) and furthermore "over"-resilience was particularly high in broadleaved forest at low elevations (interaction "ele:com" in Figure 7).
The similarity of results between the single-variable regressions and the more complex combined analysis demonstrates that interdependencies between the drivers were relatively weak and thus the singlevariable regressions were little affected by confounding with the other drivers as potentially correlated "third variables." Furthermore, the similarity of results also demonstrates that the conversion of environmental variables into 10-or 3-interval variables did not mask the major trends these drivers had on resistance, recovery, and resilience. Finally, using mean NDWI changes per pixel instead of percentage pixels with NDWI changes ≥10% yielded qualitatively the same results (Table S2).

| Resistance, recovery, and resilience of Swiss forests to the 2018 drought year
Recent studies indicate a substantial negative impact of the 2018 drought year on Swiss forests (Baltensweiler et al., 2020; The variance explained by the mixed-effects model combining the influence of different environmental drivers and their two-way interactions on the proportion of forest pixels with a decrease in NDWI ≥10% during the drought (2017-2018, resistance), after the drought (2018-2019, "negative" recovery), and the full two-year observation period (2017-2019, "under"-resilience). The bars from top to bottom in each panel are the contributions to multiple R 2 values from biogeographic regions within Switzerland to elevation, aspect, topographic position, distance from forest edge, forest composition, and interactions of these (but not region). Below region, each bar is for the three-level environmental driver or 9-combination interaction fitted with two and four degrees of freedom, respectively (see Table S1, upper part). The shading within bars is used to indicate the contributions of 1-degree of freedom contrasts: ele = contrast between lowland versus foothills and Alps, asp = north versus south, top = slopes versus valley-bottoms and ridge-tops, dis = forest edge versus non-edge, and com = broadleaf trees versus mixed and coniferous trees. The strongest shading is for the contrasts or their interactions and the least shading for the deviation from the contrasts . Additional stressors thus can diminish the natural recovery potential of forests (Dobbertin et al., 2007;Lévesque et al., 2013;Miralles et al., 2019).
For the validation of our threshold-based approach (i.e., ≥10% NDWI change), we compared our results with visual inspections of high-resolution satellite imagery (Figure 2), resulting in an overall accuracy over 83%. The producer's accuracy for forests not affected by drought was almost 99%, whereas the producer's accuracy for damaged forest was 18%. This indicates that the applied threshold only assessed severe forest damages; to include less severe damages, the threshold could be lowered. The user's accuracy was 84% and 77% for intact forest and damaged forest, respectively.
Depending on the aim of further analysis the applied threshold could be adapted. However, when we used other cutoff values of ≥5% or ≥15% NDWI change ( Figures S6 and S7 respectively) or even mean NDWI changes directly ( Figure S5 and Table S2), we obtained similar results as with the ≥10% cutoff.

F I G U R E 7
The variance explained by the mixed-effects model combining the influence of different environmental drivers and their twoway interactions on the proportion of forest pixels with an increase in NDWI ≥10% during the drought (2017-2018, "over"-resistance), after the drought (2018-2019, recovery), and the full two-year observation period (2017-2019, "over"-resilience). The bars from top to bottom in each panel are the contributions to multiple R 2 values from biogeographic regions within Switzerland to elevation, aspect, topographic position, distance from forest edge, forest composition, and interactions of these (but not region). Below region, each bar is for the three-level environmental driver or nine-combination interaction fitted with two and four degrees of freedom, respectively (see Table S1, lower part). The shading within bars is used to indicate the contributions of 1-degree of freedom contrasts: ele = contrast between lowland versus foothills and Alps, asp = north versus south, top = slopes versus valley-bottoms and ridge-tops, dis = forest edge versus non-edge, and com = broadleaf trees versus mixed and coniferous trees. The strongest shading is for the contrasts or their interactions and the least shading for the deviation from the contrasts .

| Effect of environmental drivers on forests in response to the 2018 drought
Our two types of analyses consistently showed that forest response to drought was strongly related to topographic site characteristics and forest stand structure, whereas climatic variables and interactions between drivers had comparatively small (yet still significant) effects. The same effects of drivers could be observed if they were stratified each into 10 or only 3 intervals. Because in the second case the drivers could be combined into a single analysis where they became more orthogonal to each other, we could show that their contributions to forest drought responses were largely independent of each other and, after correcting for differences between regions (used as random-effects term in our mixed-model analyses), they jointly explained a large proportion of variance in forest drought response, reducing unexplained variation to only around 20% (see Figures 6 and 7). Thus, given that the results of the combined analysis were consistent with those of the separate analyses and given the weak interaction effects in the combined analysis, we focus the remaining discussion on the results from the separate analyses for each driver.
During the drought event, less resistant forests (i.e., higher proportions of strong NDWI decrease from 2017 to 2018) were more often observed on south-than on north-facing slopes, whereas this changed during recovery from 2018 to 2019 and over the whole period from 2017 to 2019, when the proportion of damaged forest was higher on north-oriented slopes. This suggests that regarding recovery and resilience, remarkably trees of all forest types were better adapted to drought on drier and warmer slopes, which has been reported previously for mixed conifer forests (Guarín & Taylor, 2005).
Furthermore, forests in topographically exposed locations (i.e., on ridges) had both lower resistance and resilience than forests in valleys, as hypothesized. As mentioned in the first paragraph of this section, these drivers were only weakly correlated with (high) elevation and had only weak interactive effects with elevation, which in contrast had rather positive effects on resistance and resilience (see Figure 4), presumably due to positive effects of temperature.
Wind and radiation affect lager parts of the forest canopy by increasing evapotranspiration when it is located in exposed terrain, in comparison with sheltered and shaded forest stands. Another explanation is the higher buffering of soil water availability in sheltered downhill positions, whereas ridge positions have lower buffering of soil water availability and thus may cause higher water stress (Hawthorne & Miniat, 2018). This buffering has several components, such as closer proximity to groundwater at less exposed positions (Western et al., 2001) and input-independent higher recession of soil moisture at exposed locations (Hawthorne & Miniat, 2018), indicating the dominance of topography over soil properties as a driver of forest responses to droughts (Yeakley et al., 1998). Also trees in ridge locations are reported to have generally more extensive root systems (Tauc et al., 2020) and may thus be genetically better preadapted to drier conditions than trees growing in other locations (de Kort et al., 2020). Nevertheless, our analysis showed that these potential pre-adaptations were apparently not sufficient to protect those trees against the water stress experienced in 2018.
At higher elevation we found smaller proportions of forested area with severe NDWI decrease during the 2018 drought confirming findings from previous studies (Lloret et al., 2007;Lobo & Maisongrande, 2006;Rita et al., 2020). Simultaneously, strong NDWI increases were prominent in the alpine regions during the drought year, which were historically energy-rather than waterlimited (Schurman et al., 2019). With higher temperature there are more frost-free days (Jolly et al., 2005), growth conditions closer to the optimum (Klesse et al., 2018), and a considerably longer growing season (Schurman et al., 2019). According to one study, these improved growing conditions resulted in increased net CO 2 uptake and increased gross primary production in 2018 compared with 2017 in some regions . Thus, the high-temperature component of the 2018 drought event might have been beneficial in humid, cold alpine areas and have outweighed the negative effects of precipitation deficits, leading to increased growth, increased canopy water content, and thus higher NDWI values (Scherrer & Körner, 2011). Positive impacts of warmer temperatures on energylimited forest ecosystems have been found often (Dawes et al., 2011;Seddon et al., 2016). However, as discussed above, the overresistance to the drought event in 2018 was often followed by a negative recovery in 2019 at these locations, suggesting a relaxation effect from the growth spurt under conditions less limiting due to increased temperatures. Furthermore, water availability in alpine areas is highly influenced by the snowpack which was not included in the present analysis. As reported in a previous study (Gazol et al., 2017), we found increasing drought resistance on steeper slopes.
Forest stands on steep slopes are perhaps more used to reduced water availability due to faster run off.
Among the drivers related to forest stand structure, forest edges were clearly more exposed to drought damage as reflected in their low resistance values to the 2018 drought. This confirmed our hypothesis and has also been observed for other drought events (Buras et al., 2018;Laurance & Williamson, 2001). Reasons for the higher vulnerability to water stress at forest edges might be the less favorable micro-climatic conditions such as higher temperatures not just at the top of canopy but also down to the soil surface (Buras et al., 2018;Chen et al., 1993) as well as lower relative humidity (Von Arx et al., 2012) that cause higher evaporation and thus accelerate soil water stress. In contrast, further inward in a forest stand the canopy provides shade that leads to reduced temperatures and lower evapotranspirative demand (Grote et al., 2016) compared with forest edges. In general, competition for water and light should be less strong at forest edges, which might explain the higher recovery rates, but this may not have compensated the other microclimatic effects during the drought. Higher post-drought recovery rates have also been observed in tropical forests (Schwartz et al., 2019).
Whereas the linear model fit between drought resistance and forest composition was nonsignificant, forest composition was the most important driver of resilience in our study. As hypothesized, forests with a high proportion of coniferous trees were more damaged by the drought than forests with another composition in the longer term from 2017-2019. This finding corroborates earlier suggestions from fieldbased studies (Arend et al., 2021;Lévesque et al., 2013;Vitasse et al., 2019) and may be explained by the shallower rooting depths of Norway spruce, while European beech has a deep root system that reaches water in deeper soil layers (Bréda et al., 2006;Brinkmann et al., 2019).
Additionally, bark beetle outbreaks affect spruce and can weaken recovery after a drought event (Biedermann et al., 2019). Low resistance is usually followed by high recovery (Gazol et al., 2017), as we observed for example for south-facing slopes and near forest edges. The net effect is then a resilient drought response.
However, resilience can also be achieved by generally small changes in response to extreme events, that is, small negative combined with small positive changes. This has been observed in particular for diverse forests (Anderegg et al., 2018;Gazol & Camarero, 2016;Jucker et al., 2014;Pardos et al., 2021), for example, through facilitation between species (Pretzsch et al., 2013). Indeed, also in our study, more heterogeneous forests also showed this type of drought response when compared with more homogeneous forests. This confirms our hypothesis of diverse forest stands being less vulnerable to drought, presumably due to niche partitioning reducing interspecific competition.
In contrast to our hypothesis, most linear fits between proportions of forest area with strong NDWI changes and climatic variables (i.e., mean temperature anomaly, total precipitation anomaly, and CWB z-score) were largely not significant (p > .05). We did observe more damaged forest during the drought under drier conditions, but the combination with temperature showed no relations. The lack of associations between forest drought responses and climatic variables in our study may be due to the fact that the drought year 2018 was so extreme that differences in "extremeness" did not matter much or perhaps to the different spatial resolution of the input data (Allen et al., 2010(Allen et al., , 2015Brun et al., 2020). The coarse spatial resolution of the weather data suppresses important local variation especially in highly varying topography (Macek et al., 2019;Potter et al., 2013;Slavich et al., 2014), which is especially the case for Switzerland where most forest is found on hillslopes rather than on flat areas.

| Trajectories of drought-induced forest damage
Although the aggregation of data into strata makes it more difficult to follow trajectories of individual pixels, we could do so in a more general way by adding to the overall analysis a conditional analysis in which we analyzed recovery in regard of the 2018 drought for different resistance scenarios in 2018 (negative, neutral or positive NDWI changes from 2017-2018). This first post-drought year is indicative for both positive and negative drought legacy effects, which have been found to affect up to three subsequent years (Anderegg et al., 2015). A large proportion of forest pixels with low resistance (see Figure 5)

| Reliability of this study
Our analysis of satellite-based data provides a consistent and groundcomparable way to assess drought-induced forest changes over large areas . Nevertheless, several limitations and assumptions must be considered when interpreting the results. Satellitebased monitoring of drought only allows to evaluate state and changes of the upper tree crown , while possible damages happening lower in the canopy or even at stem level are not visible.
Nevertheless, defoliation leads to translucency of the canopy crown and reflectance from the understory adds (unwanted) variation to the spectral signal (Barry et al., 2008). Additionally, although the used 10 × 10 m pixel size comes close to the diameter of individual canopy trees, pixels seldom represent individual trees measured in field surveys (Vicente-Serrano et al., 2013;Zheng et al., 2021).
Clouds and cloud shadows affect measured radiances and sometimes hinder a robust retrieval of surface information. While clouds can be well masked, the identification of cloud shadows is more complicated and the separation between cloud and topographic shadow is not always possible. Our cloud-masking and gap-filling approach was conservative, resulting in a consistent exclusion of cloud shade at the cost of the additional exclusion of a small proportion of forested areas shaded by topography.
The used NDWI is sensitive to canopy water content, as demonstrated by our sensitivity analysis. However, observed NDWI changes do not only represent drought effects on canopy water content but may also reflect dynamics in canopy leaf area (Gao, 1996). It must be noted that leaf area can also be a result of land cover change (harvest) and damages caused by storms, pathogens, and bark beetle (Ips typographus L.) outbreaks (Hotovy & Jenicek, 2020;Olthof et al., 2003;Wood et al., 2018). Because we did not disentangle the contribution of these individual factors on observed NDWI changes, we cannot claim causality between drought and forest damages in all cases.
For our study, we took 2017 as a reference year, although summer temperatures were already high compared with a longer reference period (MeteoSchweiz, 2017). The contribution of possible pre-stress could have increased vulnerability and caused legacy effects from earlier events (e.g., the drought in 2015) for which trees might on average take 2 to 4 years to recover (Anderegg et al., 2015).
Such effects cannot be separated with only one base year. A longer time series including different sensor systems such as Landsat imagery at 30 m spatial resolution would be needed to increase the reference period and account for potential legacy effects. Furthermore, we only looked at the first year after drought and thus could not assess potential longer-term recovery or decline, including additional stressful or beneficial effects during later years. There is a possibility that not only growth responses but also canopy water status recovery might have continued beyond 2019. If the time span between drought events becomes shorter than the recovery phase of forest ecosystems, permanent damage is to be expected (Schwalm et al., 2017). We plan to apply the methods introduced here for the analysis of further years to potentially detect long-term legacy effects of the 2018 drought event.
We applied a 20 m inward buffer on the forest mask to avoid mixed pixels along the forest border. This filtering probably resulted in an underestimation of the total forested area affected by drought because drought tends to cause greater damage at forest edges. The applied forest mask is from 2020 and land-use changes in forest areas are thus included in the analysis. The only solution to this problem would be to generate a forest mask from the spectral imagery for each of the years 2017, 2018, and 2019 and restrict the analysis to the areas classified as forest in all years.
Additionally, the forest mask does include smaller clearings (e.g., in the north-east and south-west of the forest area in Figure 2d,e,f) and forest trails resulting in mixed signals that biased the NDWI maps. Particularly the drought response of grassland as depicted by NDWI changes is more extreme than the response from forests; and mixed spectral signals result in reduced resistance and resilience values and higher recovery values. Furthermore, mixed signals due to understory could not be erased. The understory might have been less affected due to partial shade and lower temperatures (Grote et al., 2016) and thus dampened the measured decrease in NDWI or led to an underestimation of drought damage in 2018, since shallow-rooted understory may have experienced water stress earlier than deep-rooted trees. Further, mixed pixels may cause an overestimation of forest recovery due to beneficial effects of resource reallocation from stressed trees to other plants. Non-congruent spatial overlap of NDWI maps with spatial environmental variables from different data sources may reduce correlations between these variables and forest drought response.
The coarser spatial resolution of precipitation and temperature data, for example, might have dampened a possibly stronger correlation with forest drought response. Thus, data at a finer spatial scale would be needed. Locally more precise data could be obtained from measurement stations, with the disadvantage of spatial nonconsistency. Alternatively, weather data could be improved by integrating a DEM at finer spatial resolution (Fang et al., 2013;Hutengs & Vohland, 2016). Finally, the list of environmental variables used in this study is not complete and thus cannot fully explain drought-induced impacts on forest ecosystems. We did, for example, not consider soil properties such as soil water holding capacity that are directly related to plant growth (Piedallu et al., 2011) and can thus also affect and partly explain forests' response to drought events.
Finally, we used an indirect analysis approach in the way of aggregating pixels, which had NDWI responses passing a threshold of 10%. We have also analyzed the direct values of NDWI change ( Figure S5 and

| CON CLUS IONS
We conclude that satellite-derived measurements of the NDWI are applicable to assess the impact of extreme events such as droughts on the structure and functioning of forest ecosystems across large contiguous areas such as entire Switzerland. We set to further elucidate drought responses in temperate forests and increase our understanding of how such ecosystems respond to this increasingly more frequent stressor. We found that forest stand characteristics such as distance to forest edge and functional group composition (deciduous broad-leaved vs. needle-bearing coniferous species) strongly affected drought response. This suggests reducing edge length per area and conifer component as potential management strategies to avoid negative effects of an increasing frequency of drought events. We also found that forests on southern slopes were more vulnerable to drought than forests on norther slopes, suggesting that drought events may push forests in already exposed conditions beyond their "adaptive" limits, possibly requiring replanting with more drought resistant genotypes or species. The effects of individual environmental variables on forest drought response were additive (only weak interactions) and valid across regional differences in Switzerland. Extrapolation beyond Switzerland would have to be carefully tested for different forest types and site characteristics to see if our approach is applicable to other temperate biomes where drought events in combination with higher temperatures increase in frequency and become more important in damage assessments.
Further research should expand the analysis of continuous time series rather than three month composites from successive years for a better understanding of the dynamics of forest response to severe water stress. This would also enable an assessment of the optimal timing of satellite-based drought surveys considering accumulated damage and evolving stress symptoms. Additionally, a better understanding of underlying mechanisms influencing drought will be needed to develop mitigation strategies against future drought events.

ACK N OWLED G M ENTS
We thank Lars Waser from the Swiss Federal Institute for Forest, Snow, and Landscape Research WSL for the provision of the forest mask and tree type map. We also thank MeteoSwiss for the provision of temperature and precipitation data and the European Space Agency for access to Copernicus Sentinel data. We further thank the reviewers for their comments which greatly improved our analysis.
Funding was provided by the University Research Priority Program in Global Change and Biodiversity. Open access funding was provided by the University of Zurich.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
JS, AD, and MJS designed the study. JS performed the data collection and analysis and wrote the manuscript together with inputs from AD, MJS, and BS.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.bk3j9 kdd8.