Resilience of Spanish forests to recent droughts and climate change

A widespread increase in forest cover is underway in northern Mediterranean forests because of land abandonment and decreased wood demand, but the resilience of these successional forests to climate change remains unresolved. Here we use 18‐year time series of canopy greenness derived from satellite imagery (NDVI) to evaluate the impacts of climate change on Spain's forests. Specifically, we analyzed how NDVI was influenced by the climatic water balance (i.e. Standardized Precipitation‐Evapotranspiration Index, SPEI), using monthly time‐series extracted from 3,100 pixels of forest, categorized into ten forest types. The forests increased in leaf area index by 0.01 per year on average (from 1.7 in 2000 to 1.9 in 2017) but there was enormous variation among years related to climatic water balance. Forest types varied in response to drought events: those dominated by drought‐avoiding species showed strong covariance between greenness and SPEI, while those dominated by drought‐tolerant species showed weak covariance. Native forests usually recovered more than 80% of greenness within the 18 months and the remainder within 5 years, but plantations of Eucalyptus were less resilient. Management to increase the resilience of forests—a key goal of forestry in the Mediterranean region—appears to have had a positive effect: canopy greenness within protected forests was more resilient to drought than within non‐protected forests. In conclusion, many of Spain's successional forests have been resilient to drought over the past 18 years, from the perspective of space. Future studies will need to combine remote sensing with field‐based analyses of physiological tolerances and mortality processes to understand how Mediterranean forests will respond to the rapid climate change predicted for this region in the coming decades.

Quantifying, understanding and enhancing the resilience of forests to climate change are major areas of interest (Nimmo, Mac Nally, Cunningham, Haslem, & Bennett, 2015;Oliver et al., 2015), especially in the Mediterranean area which is undergoing a warming that exceeds the global trend at present and in most projections (Guiot & Cramer, 2016a).
In Spain, enhancing forest resilience to drought has become a major goal of protected area management because the forests provide valuable ecosystem services including water regulation, timber and meat provision, regulation of climate and air quality, erosion control, as well as recreational and spiritual enjoyment (de la Guerra et al., 2017;Guiot & Cramer, 2016a;Scarascia-Mugnozza, Oswald, Piussi, & Radoglou, 2000;UNEP/MAP, 2016). Strategies for climate change adaptation used and proposed include favouring mixed species stands, reducing tree density, eliminating species with high water demands and introducing resilient species (de la Guerra et al., 2017). Given, however, how recent these management plans are, no studies have evaluated their effectiveness on a national scale, with the exception of Domingo et al. (2008), who compared greening trends inside and outside protected areas. Instead, researchers have focused on other ecosystem functions and services. A study evaluating multifunctionality in European protected areas found that in Spain those areas were associated with lower timber production and climate regulation functionality (van der Plas et al., 2018).
Another study found that land use change in Europe's network of protected areas (i.e. Natura 2000), towards non-natural cover was greater than in non-protected areas, calling their effectiveness into question (Rodríguez-Rodríguez & Martínez-Vega, 2018). National-scale evaluation of the effectiveness of managing protected area for climate resilience was however absent and is therefore needed.
Multispectral imagery collected by earth observation satellites provide information on the greenness of pixels, which is linked to total canopy cover, leaf biomass and photosynthetic activity of the forest stand (Asrar, Myneni, & Kanemasu, 1989; Baret & Guyot, 1991;Carlson & Ripley, 1997;Cihlar, St.-Laurent, & Dyer, 1991;Tucker & Sellers, 1986;Zhang et al., 2003). Time-series of remote sensing data have been used to detect forest mortality (Coops, Johnson, Wulder, & White, 2006;Fraser & Latifovic, 2005;Garrity et al., 2013;Hart & Veblen, 2015;Ogaya, Barbeta, Başnou, & Peñuelas, 2015) and also canopy-level responses to drought event. For instance, Gazol et al. (2018) estimated the responses of 11 tree species to four drought events in Spain, in terms of satellite-derived canopy greenness (normalized difference vegetation index [NDVI]) and stem growth derived from dendrochronology. They recorded loss of greenness and growth during the drought events (i.e. sensitivity) and recovery following rain. The study found that recovery-to-sensitivity ratios varied greatly among the 11 species and along climate gradients: conifer-dominated woodlands in semi-arid regions were reported to be most sensitive to drought but recovered quickly while broadleaf-dominated woodlands in humid temperate regions were least sensitive to drought and recovered slowly. However, this study did not attempt to quantify longterm impacts, including the legacy of previous droughts (Anderegg et al., 2015;Gazol et al., 2018;Peltier, Fell, & Ogle, 2016).
Here, we evaluate one component of forest resilience to climate change-the ability of canopies to resist loss of leaves and/or recover leaves quickly-by looking at time-series of greenness alongside time-series of water availability. The analysis of remote sensing data over time is now possible over large scales thanks to advances in cloud computing technology and offers the opportunity to evaluate forest responses to drought in new detail. We analyse randomly sampled Spanish forest dominated by 10 species groups (i.e. forest types), inside and outside protected areas, to investigate canopy responses to the strong drought events and also the more subtle accumulation of drought stress through time. We address the following questions: Finally, we tested whether long-term trends in NDVI observed from space related to basal area trends measured in field inventories. A previous study from Spain found close agreement between ground-and remotely sensed estimates of NDVI when field estimates were made at the same spatial resolution as the remotely sensed pixel (Ogaya et al., 2015), and province-level averages of field and remotely sensed estimates of greenness were also in close agreement (González-Alonso et al., 2006). Here we extend this approach by comparing satellite imagery with field data collected by the Spanish National Forest Inventory.
Given the above-mentioned increases in aridity during the past decade and the evidence of forest die-back as well as greening in Spain, we expected to observe forest growth in some areas because of recent land abandonment and decline in areas that have recently become drier. In parallel, we expected that species-dominating areas that have always been dry, that is, south-east Spain, would be more resilient to drought than historically wetter areas. Concerning forest species' responses to drought, we envisaged that conifer canopies, specifically pines, would be less resilient to drought compared to broadleaf canopies because of reported succession advancement (Carnicer, Barbeta, Sperlich, Coll, & Peñuelas, 2013); and that non-native fast growing eucalypt would also be less resilient given their high water needs (Queirós et al., 2020). Finally, given that protected areas are supposed to be managed for resilience, we expected them to respond differently to unprotected areas while recognizing that only 20%-50% are reported to have effective management (de la Guerra et al., 2017).

| Theory and definitions of resilience
Resilience to climate change has been defined in numerous ways and measured using many different approaches, making it hard to critically compare studies (e.g. Dakos, Carpenter, van Nes, & Scheffer, 2015;Hodgson, McDonald, & Hosken, 2015;Nimmo et al., 2015;Slette et al., 2019). Here we define resilience broadly as 'the capacity of forest canopy to return to a state not qualitatively different from its pre-drought state by resisting and/or recovering' (Folke et al., 2010;Hodgson et al., 2015). More specifically, we use the NDVI and leaf area index (LAI) as greenness and leaf area indices and the standardized precipitation evapotranspiration index (SPEI) as water availability indicator . The simplest resilience concept is illustrated in Figure 1e-here we see vegetation's response to specific drought events-'resistance' is the capacity to withstand or tolerate drought, its opposite 'sensitivity' is computed here as the absolute loss in greenness during a drought event; 'recovery' and 'adaptability' are both computed here as the absolute gain in greenness following a drought event. We use a segmentation approach to detect these drought-related events (detailed in Section 2.6). We also develop an approach to quantify short-term and long-term resilience to drought  based on measuring the correlation between greenness and water availability (Figure 1a-d; detailed in Section 2.5): here, short-term resilience is assessed by quantifying covariance (Figure 1c) and long-term resilience by comparing the observed greening trend with the potential greening trend (i.e. that predicted in the absence of climate change, Figure 1d). Remote sensing approaches are unable (currently) to measure the recruitment of new species and the loss of others, which may influence resilience via adaptation processes, these recruitment processes and shift in species distributions are therefore not addressed in this work (Folke et al., 2010).

| Site selection
Protected area maps were obtained from the Spanish Ministry of Agriculture, Food and Environment. We focussed on sampling locations categorized as Sites of Community Importance (known in Spain as Lugares de Importancia Comunitaria or LIC), and Special Conservation Zones (known in Spain as Zonas Especiales de Conservación or ZEC) within the Natura 2000 network, contrasting these sites with those without legal protection (Figure 2a). These sites include private and public lands that are meant to be managed in an ecologically and economically sustainable way, including agricultural activities, hunting and tourism (Martínez-Fernández, Ruiz-Benito, & Zavala, 2015). ZEC are supposed to be well-managed multifunctional woodlands that provide different ecosystem services, however as mentioned in the introduction, ecologically functional goals are not always achieved because of lack of proper management and funds.
In 2017, 20%-50% of Natura 200 sites were LIC whose management plans were still under consideration for approval to become ZEC. Note that Nationally Designated Protected areas-which are concentrated in Catalonia and Andalusia-were omitted from the analyses because these are more strictly managed, usually prohibiting all human activity. It would be interesting to differentiate between Natura 2000 and Nationally Designated Protected areas in a future study.
Spain has the largest forest cover of any country in the Mediterranean with more than a quarter of its territory dedicated to nature conservation, and over 1,500 protected natural areas (FAO & Plan Bleu, 2013; de la Guerra et al., 2017). Stratified random F I G U R E 1 Conceptual diagram explaining studied resilience measures: (a) A regression line (purple) is fitted to the standardized precipitation evapotranspiration index (SPEI X ) time-series (darkred) representing long-term climate changes; (b) a regression line (blue) is fitted to the deseasonalized time-series of normalized difference vegetation index (NDVI; dashed green) representing long-term change in greenness; (c) covariation is estimated between detrended NDVI (green) and detrended SPEI X (red) from a linear model between the two; (d) deseasonalized time-series of potential greening (dotted green) and its corresponding trend (orange) are determined and the difference between potential and observed greening (shaded in red) corresponding to the long-term climate change influence on greenness; (e) to gain an alternative perspective, deseasonalized NDVI is segmented (dashed purple line) and two biggest negative changes representing forest canopy sensitivity (i.e. NDVI loss) are detected followed by the recovery measures (i.e. NDVI gain; solid purple segments) sampling was used to select 4,500 1-km 2 pixels representing 10 forest types and three protection statuses across Spain. In all, 10 dominant species groups (forest types) were considered for the analysis identified from European forest species distribution maps ( Figure 2; Brus et al., 2012). We sampled 450 pixels randomly per forest type from areas that were over 60% forested selecting 10 of the most abundant species groups ( Figure S2) as computed from species distribution maps by Brus et al. (2012). The dominant forest species dataset used includes broad taxonomic categories, some of which span several bio-climatic regions, such as oaks and maritime pines ( Figure 2b).
Previous land cover could influence plant-water relations within sites because former farmland is likely to be on deeper soils, and that, in turn, could influence the resilience of forests. To test its influence, the land-cover type of each of our MODIS pixels was obtained by extracting data from the 1990 and 2000 CORINE Land Cover rasters (CLC2000_CLC1990_V2018_20 and CLC2006_CLC2000_V2018_20), and then downscaling from 100 to 500 m resolution using a nearest neighbourhood method ( Figure S3: Corine 1990 and 2000 classifications for our pixels). Out of 4,500 pixels, 342 were classified as agricultural, irrigated or non-vegetated lands in 1990 and 2000 years and were excluded because of inconsistencies with the map of Brus et al. (2012). Of the remainder, 304 (8.6%) were classified as 'previously agricultural or irrigated land' in 1990 and had transitioned to forest by 2000. We included previous land cover (i.e. agricultural vs forested in 1990) as a factor in our models.  Brus et al. (2012). Note that the 'other broadleaves' category represents planted woodlands dominated by species such as almond, carob, elm, lemon, poplar, whereas 'other conifers' represent species such as cedars and junipers with tree-ring derived resilience measures . The dataset comprises multispectral reflectance measurements recorded by MODIS Aqua and Terra satellites at a spatial resolution of 250 m (Didan, 2015a(Didan, , 2015b, from which 16-day composites were constructed. These data have been corrected for atmospheric and bidirectional surface reflectance effects, and water, clouds, heavy aerosols and cloud shadows have been masked out by Nasa. We further masked time-series of NDVI following quality assurance information and aggregated to monthly values at the resolution of 500 m. LAI was also obtained from MODIS product MODIS/006/MCD15A3H

| Datasets
for the year 2013 at 500 m resolution (Myneni, Kmynnyazikhin, & Park, 2015) to relate losses and gains in greenness to biomass on the ground. Summary data were exported for analysis in R. It is known that NDVI and LAI are nonlinearly related, with NDVI saturating at high LAI ( Figure S4), and this nonlinearity is considered when interpreting results (le Maire et al., 2006;Turner, Cohen, Kennedy, Fassnacht, & Briggs, 1999). A nonlinear model was fitted to average NDVI for the year 2013 with average LAI, and the model was used to convert NDVI time-series into LAI time-series and interpret all results obtain from NDVI analysis in LAI terms as well (Table S1, M5; Figure S4).
An absolute measure of yearly climatic water balance was computed for each site to evaluate the effect of the latitudinal climatic gradient on forest response to drought. It corresponds to monthly precipitation minus potential evapotranspiration averaged across the 18 years.
Potential evapotranspiration was calculated using the 'Penman' function in the SPEI package , by inputting minimum and maximum temperature, mean daily external radiation, wind speed and pixel latitude. All climatic variables were extracted from Global Land Data Assimilation System 2.1 products (Rodell et al., 2004). In addition, time-series of monthly absolute water balance were also used to be compared with the relative water balance index, SPEI (used in this paper to evaluate the intensity of drought).
Time-series of SPEI, representing relative water availability, were extracted from a database available at 1.1 km resolution for the entirety of Spain from 1961 to 2017, and used as a climate change index . It is important to note that SPEI provides an index of temporal change in water availability within a site. It should not be used to make direct comparisons of absolute water availability among sites but can be used to compare the relative intensity of a dry or wet period. For example, in this paper, we are mostly interested in the severity of drought events relative to the severity of drought events experienced over time in a site and use this relative severity (i.e. SPEI) when comparing forest responses. SPEI is widely used as an alternative to SPI, which is based solely on precipitation (Azorin- Bottero et al., 2017;Gazol et al., 2018;Tejedor, de Luis, Cuadrat, Esper, & Saz, 2016; Vicente-Serrano, Beguería, & López-Moreno, 2010). SPEI is computed as follows: first, a monthly water balance is calculated for every month from 1961 to 2017, by subtracting weekly potential evapotranspiration from weekly precipitation. Second, because the impacts of drought can be cumulative, a sequence of alternative SPEI time-series is created from SPEI values averaged over a different number (X) of months. In this study, SPEI is calculated for X ranging from 1 to 48 months. The cumulative water balance is calculated from the weighted average of a particular month and the X − 1 previous months using a rectangular weighing function.
The SPEI X is then obtained from normalizing water balance series into a three-parameter log-logistic distribution using non-biased probabilistic weighted moments to calculate the parameters. The parameter estimation of the log-logistic probability distribution and detailed calculation procedure for SPEI index can be found in Vicente Datasets mentioned in this section, with the exception of SPEI, were all available through Google Earth Engine.

| Greening and climate change trends
A particularly effective method of monitoring vegetation from space is detecting and understanding changes in time-series of greenness (Jamali, Jönsson, Eklundh, Ardö, & Seaquist, 2015). These changes tend to fall under three categories: seasonal changes, gradual changes and abrupt changes (Verbesselt, Hyndman, Newnham, & Culvenor, 2010).
Recent methods of analysing these changes consist of decomposing the time-series into 'trend' which captures the gradual changes caused by interannual climate and land use variations as well as the abrupt changes cause by disturbances, 'seasonality' which follows annual temperature and rainfall, and 'remainder' variations unexplained by the former factors (Jamali et al., 2015;Verbesselt et al., 2010). In this work, we focus only on gradual changes in NDVI and LAI trend, which are the type of change that drought usually causes. In this first part we deseasonalized the NDVI time-series to relate it to non-seasonal climate change index SPEI using the DBEST package in R (Jamali et al., 2015;Jamali & Tomov, 2017). The DBEST function detected abrupt changes in 20 of the 4,158 pixels, which were excluded from subsequent analyses as they represent extreme events such as fire, clear-felling and sensor-related artefacts which are beyond the scope of our study. A further 952 plots were eliminated because of missing NDVI or SPEI data-points leaving 3,528 plots for further analysis. We used two linear models (Table S1,  (1) are normally distributed residuals. Subscripts in these models indicate C for climate and G for observed greenness. is the predicted initial SPEI X or NDVI of a plot (i.e. on February 2000), and represents the long-term linear trend in climate (i.e. SPEI) or greenness of a plot, respectively.
Further linear models were used to evaluate how greening trends vary across Spain, testing whether they vary with elevation, average water balance, average greenness, forest type and protection status (Table S1, model M10-M11). Models were conducted using spatial simultaneous autoregressive error estimation (SSAEM), which account for the residual spatial structure within the neighbourhood of pixels (spatialreg package in R; Bivand & Piras, 2019). AIC and correlograms at varying neighbourhood distances were used to determine the best-supported model to account for spatial autocorrelation.

| Determining the drought accumulation period across Spanish forests
To determine the drought accumulation period over which drought affects forest greenness, we recorded which SPEI X time-series correlated most closely with the forest greenness time-series at a site.
We used a cross-correlation function with zero time lag (as the X scales of SPEI are already computed to account for time lags) correlating detrended SPEI X against detrended, deseasonalized NDVI for drought accumulation periods ranging from X = 1 to 48 months, and the value of X that gave rise to the highest correlation was noted: we refer to this as the drought accumulation period, and this period was used in all subsequent analyses. This allowed us to account for lags in forest response to drought and delayed mortality to which some species are susceptible, and for different water reservoirs that cause a delay in water deficit and surplus effects (Vicente-Serrano et al., 2013).

| Time-series analysis of forest response to drought
SPEI X and NDVI time-series were used to quantify two metrics for each pixel: short-term covariance and long-term climate change influence (CCI; see Figure 1). A linear model (Table S1, M4) was fitted to the deseasonalized NDVI time-series using smoothed SPEI X time-series, where deseasonalization was achieved using the DBEST package, and smoothing using locally weighted smoothing (θ = 0.1) in R, to eliminate noise and make SPEI X comparable to the NDVI trend: where t is time (1-215 months) within the 18-year time-series, G , G and G are coefficients estimated by least squares regression and G is the residual. Subscripts in these models indicate P for potential greenness and G for observed greenness.

Substituting Equation (1) into Equation (3) we get:
This is reorganized to give: Here P is the potential greening trend; it is comprised of the observed greening trend P and the influences of long-term changes in climate ( G . C ) on greenness. G provides an estimate of how strong the covariance between greenness and changes in water availability (CGW) is G − P = G . C corresponds to estimated long-term CCI on greenness.
G is the residual which contains the nonlinear trend, the seasonality and the noise (also normally distributed).
Regression modelling produces two new summary variables for each pixel: short-term CGW and long-term CCI. It is key to understand that CGW, represented by the coefficient 'γ' in Figure 1c which quantifies how greenness varies with drought index on a monthly basis, relates to both losses and gains and is not an indication of resilience alone. Large CGW values correspond with big losses and/ or big gains. Strong and weak covariance could indicate different mechanisms to overcome droughts better known as drought avoidance and drought tolerance mechanisms, respectively, which relate to anisohydric and isohydric behaviours (Hwang et al., 2017;Roman et al., 2015;Sade, Gebremedhin, & Moshelion, 2012). CCI, on the other hand, represents the long-term CCI on greenness; it can be represented as the difference between the actual greenness and the potential greenness as demonstrated in Figure 1d if greenness in the long term is affected to similar degree to the short term by the coefficient 'γ'.
Statistical modelling was used to evaluate how short-term covariance and long-term CCI varied across Spain, in relation to elevation, average water balance and average greenness (Table S1, M14-M21).
Additional models included forest type and protection status as explanatory variables.
We evaluated whether field data confirmed remote sensing analyses at the national scale using the Spanish National Forest Inventory, albeit with plots that differed in resolution from the satellite data. The Spanish Forest Inventory is a systematic network of plots each km 2 of forested areas in Spain (Villaescusa & Díaz, 1998;Villanueva, 2004). For basal area calculations, we used the permanent plots between the third (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) and fourth (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) (4) = P . t + P + P (t) , restricted our analyses to pixels in which forest was relatively homogeneous. To do this, a cloud-free image of Spain was obtained (by mosaicking Sentinel-2 level 1C images from May to August 2016) and NDVI computed for pixels of 25 and 500 m resolution centred on the plot locations. Sites which had a difference in NDVI over 0.01 between the two resolutions were excluded from further analyses leaving 2,150 out of 4,649 pixels/plots to analyse.
Modis NDVI data were extracted using the approach described in Section 2.3.1 and the absolute NDVI trend (determined from the DBEST package) was computed between the dates of the third and fourth national forest inventory. Spatial autoregressive modelling was used to determine the relationship between the greenness and basal area change.

| Breakpoint analysis approach to resilience estimation
Breakpoint analyses were conducted to detect the greatest changes in NDVI and LAI within each time-series, from which loss and gain indices were estimated. The time-series NDVI and LAI were decomposed into trend, seasonality and remainder (see Section 2.3.2), and the trend segmented using the DBEST package in R (Jamali & Tomov, 2017). As illustrated in Figure 1e, loss and gain were estimated from the two most negative segments followed by a positive segment, respectively.
Loss was estimated as the mean of the absolute negative changes in NDVI (i.e. NDVI peak − NDVI pre-event ) and gain as the mean increase in NDVI following these negative changes (i.e. NDVI post-event − NDVI peak ).
Minimum water availability (minimum SPEI) as well as maximum water availability (maximum SPEI) were extracted over those two detected loss and gain periods, respectively, and their mean computed for each series, representing water deficit and water surplus events, respec- To determine the relationships between loss and gain, and to test whether it varied systematically among forest types and protection statuses, we used ordinary least-square regression. Loss and gain were log-transformed to improve residual normality (Table S1, model M5; Figure 10a; Figure S6). Predictions were then back-transformed to arithmetic scale with a correction factor (Figure 10a; Figure S6; Baskerville, 1972).
To help elucidate the drivers of resilience, losses and gains in greenness and leaf area, and their ratio estimate, were modelled as a function of elevation, average water balance, average greenness, relative water deficit (-minimum SPEI), water surplus (maximum SPEI), forest type and protection status again using SSAEM (Table S1, M22-M39).

| Water balance is the biggest determinant of forest greenness and leaf area
As expected, forests were greenest in the wetter regions of Spain

| Long-term trends in water availability and greenness
Over the past 18 years, 73% of Spanish land has become drier, with an average change in SPEI 1 of −0.24, indicating the mean, which is almost a quarter of a standard deviation drier than the average for 1961 to 2014 ( Figure 4). This corresponds to an increase in water deficit of about 2 mm/year as calculated from the linear relationship between SPEI 1 and water balance (Figure 4c). Of the forest pixels that we randomly sampled over Spain (Figure 2b), 35% showed a significant drying trend (linear modelling, p < .05), 37% did not change significantly and 27% got wetter (Figure 5a). The climate of the wetter regions changed least while the drier regions became increasingly arid (i.e. Figure 4b,c).
Despite these strengthening water deficits, 75% of Spanish forest pixels became greener over time and only 12% became less green (Figure 5b). On average there was a 11% increase in LAI (from 1.7 to 1.9 m 2 /m 2 ), and some variation with canopy density (i.e. average LAI): less green sites (lower quartile LAI = 1.36) accumulated 0.1 m 2 /m 2 of leaf area while the more green sites (upper quartile LAI = 2.26) accumulated 0.3 m 2 /m 2 over 18 years ( Figure 6).
Correlation of remote-and plot-based observations of forest change (see Section 2.5 for details) showed a significant positive logarithmic relationship between basal area and average deseasonalized NDVI at the time of the third inventory, with a Nagelkerke R 2 of .22 (Figure 3a). A significant positive linear relationship was found between basal area change and NDVI trend change between the two inventories ( Figure 3b) with a Nagelkerke R 2 of .12. These analyses build confidence that our remote sensing analyses of resilience correspond with changes observed on the ground.
Examining individual time-series, we found a close covariance of greenness with climatic conditions: SPEI X explained 49% of the variation in canopy greenness (average R 2 of 3,182 linear models fitted to paired SPEI-NDVI time-series). Only 19% of the pixels responded immediately to water deficit (i.e. correlated best with a drought accumulation period X of only 1 month).
In the long term, CCI on greenness was significant in 58% of the plots (with 56% negative effects and 2% positive effects, p < .05) but changes were small in magnitude (Figure 7). The CCI on NDVI varied across Spain as follows (Figure 7a CCI [LAI] increased from −0.020 at the driest sites to 0.017 at the wettest sites (i.e. from −1,181.5 to −301.5 mm, these representing the 5th and 95th percentiles of observed distribution in annual water balance). It also increased with LAI from −0.013 at 0.88 to 0.011 at 2.93 (again the 5th and 95th percentiles of LAI). Again, this means that the CCI on greenness caused the LAI trend to be significantly smaller in drier or more dense sites (Figure 7b).

| Water availability is the main driver of shortterm changes in forest greenness
The short-term effect of climatic variation on canopy greenness (CGW, see Section 2.5) varied across Spain as follows (Figure 8) (coefficients in bold are significant at p < .01): The terms in square brackets are scaled environmental variables.
The model demonstrates that the covariance was weaker in elevated areas. The covariance was also weaker in forest with high NDVI but  and gains in LAI were greatest in dry regions (Figure 9(1)). The relative water deficit during a drought event (-minimum SPEI) and relative water surplus post drought (maximum SPEI) both had a significant effect on all resilience components. The canopy greenness losses and gains were greater with an increase in water deficit ( Figure 9a,b(2)). Furthermore, pixels hit by stronger drought were more resilient as deduced from the water deficit negative relationship to the gain-to-loss ratio (Figure 9c(2)). The NDVI gain, as well as gain-to-loss ratio, were equally affected by how wet the period that followed the drought was (Figure 9b,c(3)). Denser forests (i.e. forest with higher LAI) lost and gained more due to droughts ( Figure 9(5)). Forests at higher elevation were more resistant to drought events (Figure 9a(5)).

| Influences of forest type on drought response
Forests exhibited a spectrum of responses to drought, with breakpoint analysis indicating that some forest types lost and re-gained little NDVI or LAI while others fluctuated greatly ( Figure 10a). Gains and losses were correlated across forest types (r = 0.66 for NDVI and r = 0.48 for LAI; Figure 10a).
Specifically, Castanea-dominated forests had below-average losses and gains indicative of drought tolerance while maritime pines had above-average losses and gains indicative of drought avoidance (Figure 10a). This range of responses was confirmed by the analyses of the entire time-series; there was a close correlation between the short-term covariance of F I G U R E 5 Changes in Standard Precipitation-Evapotranspiration Index (SPEI) and greenness over the past 18 years in studied Spanish forests: (a) Relative frequency histogram of SPEI X change over forest plots with significant changes in red and non-significant changes in grey (linear modelling, p < .05); (b) relative frequency histogram of changes in greenness over forest plots with significant changes in pink and non-significant changes in grey (linear modelling, p < .05). Black lines indicate average of change in SPEI and normalized difference vegetation index F I G U R E 6 Relationship between plot-based measurements of basal area and remotely sensed measurements of greenness. (a) Basal area measured during the third Spanish national forest inventory as a function of average deseasonalized normalized difference vegetation index (as extracted from DBEST) for the same date; (b) greenness change as a function of basal area change between the third and fourth Spanish national forest inventory greenness and SPEI, and the greenness gains and losses determined by breakpoint analysis within the different forest types ( Figures S13 and S14). Species groups that lost and gained little after drought had a lower short-term CGW while species groups that lost and gain a lot had higher short-term CGW (indicated by circle sizes in Figure 10a); short-term covariance F I G U R E 7 Climate change influence (CCI) on greenness (i.e. difference in greenness trend cause by long-term trends in climate), as a function of elevation (m), water balance (mm/year) and average LAI (m 2 m −2 year −1 ). (a) CCI as computed from modelling normalized difference vegetation index with Standard Precipitation-Evapotranspiration Index (SPEI) as function of environmental variables; (b) CCI as computed from modelling LAI with SPEI as function of environmental variables. Uncertainty around regression lines is shaded in the environmental variable respective colour F I G U R E 8 Variation of short-term covariance between greenness and Standard Precipitation-Evapotranspiration Index (SPEI) along forest density and elevation gradients. (a) Covariance between normalized difference vegetation index and SPEI; (b) covariance between LAI and SPEI. Uncertainty around regression lines is shaded in the environmental variable respective colour F I G U R E 9 Effects of environmental conditions on (a) loss during drought; (b) gain after drought and (c) gain-to-loss ratio of LAI. The effect sizes shown are coefficients (mean ± 99% confidence intervals) estimated by spatial autoregressive modelling. Explanatory variables are scaled, so coefficients indicate the effect on these components of shifting to an environment that is +1 SD from the mean environment. The dashed lines indicate the average loss, gain and gain/loss ratio for the mean environment is correlated to loss (r = 0.87 for NDVI and r = 0.48 for LAI) and gain (r = 0.44 and r = 0.32 for LAI, p < .05 in both cases) ( Figure 10a).
We found that the non-native Eucalyptus had low gain-to-loss ratios following droughts (Figure 9b), which explains the decline in leaf area in 26% of examined Eucalyptus plots, a significantly higher percentage than the average of 12% (Table S2). Looking at the long-term CCI, we found that growth in most Spanish forests is similarly impacted by SPEI trends (Figure 10c). Only two groups, Pinus pinaster associated with xeric sites, as well as Quercus robur and Quercus petraea, were more severally pressured by climate change (Figure 10c). Combining this with previous findings, we conclude that the overall greening of Spain is possible in part because of species which are resilient to drought while the decline in many Eucalyptus plantations (26%) is mostly caused by their low short-term resilience.
F I G U R E 1 0 Summary coefficients for the forest types of the spatial autoregressive model performed to evaluate resilience components. (a) Relationship between the loss, gain and the short-term response to drought of forest types; (b) gain-to-loss ratio (log scale) and (c) longterm climate change influence on greenness. Value for forest type indicates mean response for forests outside of protected areas. Overall mean for forest outside protected areas indicated by the black dot in (a) and the dashed lines in (b) and (c), and the error bars represent 99% confidence intervals. Note losses and gains were log transformed to improve normality of residuals, since both measurements are right skewed, and because gains tend to be smaller when losses are large F I G U R E 11 Leaf area losses and gain in the short term across different protection statuses. (a) LAI losses and gains to sever droughts with fitted line back transformed from the log-log relationship between the two; (b) summary coefficients for the protected areas of the spatial autoregressive model performed to evaluate short-term covariance between Standard Precipitation-Evapotranspiration Index (SPEI) and normalized difference vegetation index and (c) summary coefficients for the protected areas of the spatial autoregressive model performed to evaluate short-term covariance between SPEI and LAI. Values in (b) and (c) are for all forest types. Averages for non-protected areas are indicated by the dashed lines and the error bars represent 99% confidence intervals 3.5 | Resistance of protected areas to climate change and effect of previous land cover Canopy greenness and leaf area inside protected areas were more stable over the 18-year time-series (Figure 11c). Looking at the covariance between SPEI and canopy greenness results, protected forests were significantly less sensitive to drought than unprotected forests (i.e. lost and gained less greenness during droughts; Figure 11b). On the other hand, the equivalent response in leaf area change was not significant indicating that the lead area did not change much in these forests (Figure 11c). No significant difference was found in the gainto-loss ratio in the most extreme droughts yet, since, as mentioned previously, Spanish forest canopy greenness was found to typically resilient (Figure 11a).
Including whether or not forested areas were previous agricultural lands did not improve the fits of any of our models (i.e. changes the AIC were less than 2; Table S1).

| Long-term and short-term CCI on forest canopy
Regression and breakpoint analyses of the NDVI time-series shed light on why Spanish forest canopies are so remarkably resilient to climate change, greening despite the region becoming drier. The NDVI time-series covaried closely with relative water availability (R 2 = .49) which agrees with previous findings (Gouveia et al., 2012).
In the short term, this covariance caused significant variation in canopy greenness and leaf area in response to drought events, with average losses in LAI during strong droughts amounting to 0.32 m 2 / m 2 ; however, in the long term, this covariance had a much smaller influence on greening trends with predicted changes to LAI amounting on average to 0.0002 m 2 m −2 year −1 with a standard deviation of 0.004. Drought effects on vegetation are commonly defined as multi-scalar phenomenon, progressing from stomatal closure, and/ or leaf abscission to embolism and hydraulic failure in response to water deficit (McKee et al, 1993 ;Vicente-Serrano et al., 2010), but few studies on forest resilience have dealt with this temporal aspect of trees response to drought Schwalm et al., 2017). The use of SPEI at different accumulation periods X allowed us to account for lags in forest response to drought and to better capture the effect of water deficit on vegetation. Recent analyses of Spanish forest indicate that drought was responsible for 50% of damage to forest canopies ('damage' defined as >25% defoliation of the canopy), followed by insects and pathogens that were responsible for 24% of damage (IDF España, 2017). However, these analyses focus on the sensitivity of woodlands not recovery. On average, 80% of losses in greenness were recovered in the 18 months following a drought event, while the remaining 20% were recovered in subsequent years to give the long-term greening trend we observe (Figure 10a,b; Figure S7). These results suggest that forests in Spain mostly undergo defoliation that is easily recoverable in the wet period following a drought while 20% of the damages that can be more permanent are outgrown in the next couple years before another drought happens. Water availability in period after severe drought is important for forest recovery (Jiang et al., 2019), as we indeed found that greater water availability (maximum SPEI) led to greater recovery while minimum SPEI was associated with bigger losses in greenness and leaf area.
Surprisingly, forests exposed to the strongest droughts had the highest gain-to-loss ratios (Figure 10b). This phenomenon was observed in several other studies on the effect of drought on forest canopy and basal area (Delissio & Primack, 2003 (Delissio & Primack, 2003;Navarro-Cerrillo et al., 2019;Schwartz et al., 2019;Sohn, Hartig, Kohler, Huss, & Bauhus, 2016). In turn, the reduction of competition aboveground and belowground can boost the recruitment of more resilient trees and the productivity of the surviving trees during the wet period following an intense drought (Schwartz et al., 2019;Slik, 2004). It is possible however that observed increases in remotely sensed short-term resilience are not related to the detected increases in basal area increment (Dorman et al., 2015;Gazol et al., 2018), but instead to the increase in understorey photosynthetic activity (Breshears et al., 2005;McDowell et al., 2008;Schwartz et al., 2019).

| Dense canopies are most sensitive to drought
Changes in satellite-detected greenness are known to be correlated to leaf biomass change and photosynthetic activity (Asrar et al., 1989;Baret & Guyot, 1991;Carlson & Ripley, 1997;Cihlar et al., 1991;Tucker & Sellers, 1986;Zhang et al., 2003), and many studies have used this information to estimate forest biomass from satellite (Galidaki et al., 2017;le Maire et al., 2011). Examining the relation between the Spanish national forest inventory basal area and greenness we found a significant positive relationship which supports these studies. Consequently, our results indicate that in Spain higher canopy density increases competition and transpiration leading to increased sensitivity to drought (i.e. losses in LAI) which agrees with previous studies on forest density effects (Bottero et al., 2017;Navarro-Cerrillo et al., 2019;Raz-Yaseef, Yakir, Schiller, & Cohen, 2012). We found that forests that have denser canopies, that is, have a high average NDVI or LAI, are more sensitive to changes in SPEI. The effect of forest density becomes however insignificant when looking at the gain-to-loss ratio to extreme drought (shortterm CGW), meaning that density-caused losses in greenness or leaf area are always followed by similar NDVI or LAI gains. We believe that the lack of forest-density effect on short-term greenness resilience is mostly due to the overall resilience of the forests, and that stronger droughts will cause this effect to be significantly negative.

| Differential response of forest types and protected areas
Spanish forests are currently dominated by pines, which occupy 67% of forested lands, including P. pinaster and P. sylvestris in big proportions ( Figure S2; Brus et al., 2012). The second most dominant species groups are oaks (22% including Q. robur and Q. petraea), secondary successional species, which are favoured by managers and are reported to be advancing on pines (Brus et al., 2012;Martín-Alcón et al., 2015;Pausas et al., 2004). Understanding the resilience to drought of these two groups is therefore paramount for Spanish ecosystems. Here, forest types differed in response to droughts, indicating differences in the drought resilience behaviours of species that dominate these forest types (Hwang et al., 2017). We found that Castanea lost little greenness during drought, consistent with anisohydric species that keep their stomata open and continue photosynthesizing for longer when droughts occur (Martínez-Sancho, Navas, Seidel, Dorado-Liñán, & Menzel, 2017;Mota et al., 2018).
This 'tolerance behaviour' results in their leaf water potential becoming increasingly negative during drought periods, which can push anisohydric trees beyond their limited hydraulic safety margins, leading susceptible species to hydraulic failure, that can be temporary or permanent (Kannenberg et al., 2019;Martínez-Sancho et al., 2017). At the other end of the spectrum, P. pinaster lost many leaves in response to drought. Isohydric species, like P. pinaster, tend to close their stomata to preserve hydraulic conductivity and then drop leaves as the drought continues (Galiano, Martínez-Vilalta, & Lloret, 2011). Species with an 'avoidance behaviour' can die during prolonged droughts as a result of hydraulic failure. Some studies also suggest that carbon starvation can occur (Galiano et al., 2011;Savi et al., 2019;Sevanto, Mcdowell, Dickman, Pangle, & Pockman, 2014), but other researchers remain unconvinced that this mechanism is important (Körner, 2013;Muller et al., 2011).
Plantations of Eucalyptus responded differently to other forest types, recovering slowly from drought ( Figure 10a) and in a number of cases (26%), higher than that of other species (12%), failing to recover (Table S2). There are reports of Eucalyptus declining in response to drought in the Mediterranean regions of their native Australia (Brouwers, Matusick, Ruthrof, Lyons, & Hardy, 2013;Matusick et al., 2018;Matusick, Ruthrof, & St Hardy, 2012). In Spain, Eucalyptus is planted in wetter regions that have not become drier over the past 18 years, and most reports and studies discuss Eucalypt in the context of timber production, pests and diseases (Queirós et al., 2020 ;IDF España, 2017) rather than their response to drought (Barradas et al., 2018;Peñuelas, Lloret, & Montoya, 2001).
The study provides no evidence of widespread permanent hydraulic failure of either isohydric or anisohydric species, as most canopies we investigated recovered in greenness following drought events. With the exception of Eucalyptus plantations, there were no significant differences in the gain-to-loss ratio of forest types, which might indicate that these contrasting behaviours were similarly effective in terms of canopy resilience. These findings need to be treated with caution, given that few of the forest types were comprised of single species. For instance, P. sylvestris is recognized as an isohydric species (Galiano et al., 2011), but forest classified as dominated by this species did not respond to drought in the way expected of an isohydric species-they lost about average canopy greenness during drought (Figure 11a)-so we suspect that other woody species contributed to the greenness signal.
Additionally, stands with medium-low canopy cover that leave exposed understorey and shrub or herbaceous which can contribute to the greenness. Focusing analyses on stands with known species composition will help clarify the responses of species and communities to drought (e.g. Hwang et al., 2017). Furthermore, studies that link canopy greenness trends with detailed physiological measurements are needed to understand the mechanistic explanations behind resilience.
Our regression analyses demonstrate that managements of woodlands in protected areas have made their canopy greenness, that is, NDVI, more resilient to drought because their short-term responsiveness to water availability is relatively low (Figure 11b).
When looking at LAI, however, the response to water availability was not significant different in protected areas, contradicting the NDVI results (Figure 11c). To extract more information on the matter, we decided to look at the number of declining LAI pixels instead of the effect on all pixels; it was significantly smaller in Natura 2000 forests and the decline was also significantly less pronounced (Table S2). These results indicate that Natura 2000 ZEC might have a more positive effect on the forest canopy resilience.
Field studies in the region have demonstrated that thinning can improve resilience (Navarro-Cerrillo et al., 2019) but thinning can also decrease leaf area. This is the first study to determine the effect of these protective measures on forest canopy resilience to climate change at the national scale from satellite imagery. More studies are needed to determine which measures in these protected areas are actually enhancing forest resilience. As mentioned in Section 2.2, the management of Natura 2000 sites varies greatly depending on the responsible administration, the funds available, the biogeographic region and whether or not the site is a Nationally Designated Protected area, therefore, future studies can take these factors into account to uncover more information on management effectiveness.
Finally, gain-to-loss ratios obtained by breakpoint analysis, as well as LAI response, could be failing to show that protected forests are more resilient to drought mainly because of the overall resilience of Spanish forest canopy greenness to past changes in water availability. Alternatively, given that the short-term NDVI responsiveness to drought was not very strong, it could be that protected areas have failed to make a difference in forest resilience thus far; more data are needed to evaluate this. As mentioned in van der Plas et al. (2018), many protected areas have been established recently; thus, there has been insufficient time for the implementation of their transformative management actions. In addition, many have been created to protect and restore degraded and abandoned lands which require extra time to recover (van der Plas et al., 2018). Although we found no effect of land-cover change from 1990 to 2000 on the resilience of forest canopies, we did not evaluate the effect of past forest management, which would also be interesting to study.

| Research and management implications
Returning to the questions posed at the beginning of this study, it is now possible to state that canopy greenness of Spanish forests have been remarkably resilient to droughts in the 2000-2017 period: these secondary forests are continuing to accumulate biomass and becoming greener despite climate change making this region hotter and drier. Unlike our expectations, most forest types, regardless of their drought resistance behaviour and management, recovered from losses incurred by extreme drought events. Finally, we demonstrate that protection and management of these forested lands has the potential to be effective in alleviating the effects of climate change on canopies by increasing the forest capacity to resist during a drought.
Breakpoint and regression approaches provide complementary insights into resilience and its components. Short-term covariance between water availability and NDVI or LAI was highly correlated with sensitivity and to a lesser extent recovery measures and could be used to assess the stability of a forest in face climate variation, as well as provide insights into the prevailing drought resistance mechanism of a species. The long-term CCI, on the other hand, can be used to estimate drought pressure on these forests. Finally, breakpoint analysis, which allows losses and recovery of canopy greenness due to drought events to be determined, could be used alongside forest die-off information to determine the tipping point of forests.
Our study has focused on the resilience of forest canopies to drought but shifts in species distribution and changing disturbance regimes will also influence resilience as the climate warms. Several studies have already analysed past shifts in species distribution as well as modelled and predicted future changes in plant communities and their distributions using ground data (García-Valdés, Zavala, Araújo, & Purves, 2013;Lines, 2012;Lloret, Martinez-Vilalta, Serra-Diaz, & Ninyerola, 2013;Rabasa et al., 2013;Ruiz-Labourdette, Nogués-Bravo, Ollero, Schmitz, & Pineda, 2012;Valladares et al., 2014). Most of these studies mentioned a decrease in the climatic suitability in the future for most species and predicted altitudinal shifts. A few mentioned an increase in habitat for broadleaves at the expense of pines while one study predicted an increase in the range of low altitude pines.
Furthermore, extreme drought events increase the risks of fires and insect outbreaks (Allen et al., 2010). Our methodology can be applied to such secondary events if data such as fire occurrence, insect distribution and physiology, and geology is available; and it would be the next natural step to identify management options that would maximize the resilience of forest to different extreme biotic and abiotic events. For instance, fire occurrence data can be tracked from satellite (Justice et al., 2002;Turco, Herrera, Tourigny, Chuvieco, & Provenzale, 2019), and pathogen damage can be detected from hyperspectral data, lidar data or both (Lin, Huang, Wang, Huang, & Liu, 2019;Stereńczak et al., 2019).
Furthermore, species composition and stand development information are becoming more and more accessible thanks to hyperspectral imagery and repeat lidar scans Jucker et al., 2018;Nunes, Davey, & Coomes, 2017;Simonson, Allen, & Coomes, 2012;Simonson, Ruiz-Benito, Valladares, & Coomes, 2016). Future work could also explore whether stand height, age, soil type and species composition influence resilience; and include field data at the same resolution to gain a better understanding of the linkages between remotely sensed greenness and forest change.

S.K. is supported by Cambridge Commonwealth, European
& International Trust. We acknowledge MAPA (Ministerio de Agricultura, Pesca y Alimentación) for the access to the Spanish National Forest Inventory. We thank Dr. Paloma Ruiz-Benito for helping us access the Spanish National Inventory data, Dr. Tom Swinfield and Dr. Javier Igea for their advice on spatial auto-correlation models, Dr. Thiago Burghi for his advice on modelling, and Miss Eleanor Hammond for reviewing the manuscript. Finally, we would like to thank the anonymous reviewers for the valuable comments and suggestions which helped us improve the manuscript.