Identification of forest vulnerability to droughts in the Iberian Peninsula

The increase in frequency, severity, and duration of droughts poses as a serious issue to the management of forests in the Iberian Peninsula, with particular emphasis on the decline of forest growth and forest dieback. Hence, the adoption of adaptation and mitigation measures in forest ecosystems that are more vulnerable to drought is a pressing matter that needs to be addressed in the near future. This work aims at identifying the regions in the Iberian Peninsula where forest exhibit high vulnerability to drought conditions. To accomplish that, a vulnerability map is produced by considering three pillar components: exposure, sensitivity, and adaptive capacity to drought. Exposure is estimated based on the multi-scalar drought index Standardized Precipitation-Evapotranspiration Index (SPEI) and aridity, while the remotely sensed Vegetation Health Index (VHI) and mean forested cover are used to assess the regions’ sensitivity to drought. Finally, elevation, water table depth, fire radiative energy, and annual solar irradiation are compiled as indicators to assess adaptive capacity. Principal component analysis was then applied to the three pillar components to identify the areas more vulnerable to drought. This approach allows for the identification of forested areas vulnerable to drought in terms of vulnerability classes automatically determined. Forests presented very high vulnerability in eastern Spain, and central Portugal. Within the most vulnerable vegetation communities, mosaic tree and shrub types revealed to be extremely vulnerable to droughts in the Iberian Peninsula, followed by needle-leaved forests (in Central Portugal, and Northeast Iberia). This work highlights the regions and primary vegetation communities to which the effort of adapting and mitigating drought consequences should be utterly enforced by the responsible authorities.


Introduction
Some forest ecosystems are exceptionally rich in terms of biodiversity, representing an essential part of the carbon cycle, and are crucial to the socio-economic development of regions that depend on the processing and trade of forest goods and services (Palahi et al. 2008). Countries such as Portugal and Spain, located in the Iberian Peninsula (IP), are characterized by different types of forests, which according to national statistics occupy about 36% of the Portuguese territory (Uva et al. 2015) and 39% of the Spanish territory (2019), and which exploration greatly contributes to both countries' net income (Rego et al. 2013;Pinho 2014;Uva et al. 2015;Sierra-Pérez et al. 2015). Indeed, business volume from the forest-based industry and silviculture represented 5.2% of the 2020 gross domestic product (GPD) in Portugal (https:// www. ine. pt). Furthermore, according to Eurostat (https:// ec. europa. eu/ euros tat/), Spain is the second country in the European Union in terms of absolute forest surface, only surpassed by Sweden, and the second country with an increase in forest area in the last 20 years-more than 34% of area, only behind Ireland.
Forests are susceptible to climate constrains and particularly to extreme weather events, such as droughts, heatwaves, and wildfires (Teuling et al. 2010;Choat et al. 2012;Gouveia et al. 2017;Boer et al. 2020;Anderegg et al. 2020). Particularly, drought episodes are a recurrent feature of climate that have become more frequent, more intense, and spanning larger periods of time in the last decades over the IP (Spinoni et al. 2015a, b;Páscoa et al. 2017Páscoa et al. , 2021. As a climate change hotspot (Giorgi 2006), this region is further expected to endure increasingly frequent, severe, and larger droughts in the future (Spinoni et al. 2018;Naumann et al. 2018). Changes in the frequency of droughts, which can be further aggravated by insect outbreaks and wildfires, may impose great losses to commercial forestry (Kirilenko and Sedjo 2007). The IP is characterized by large spatial differences in drought features, with more frequent episodes in the northern regions when compared with centre and south IP, but more intense and long-lasting in the southern part of the Peninsula (Domínguez-Castro et al. 2019). The influence of droughts in different species of pines and locations has been extensively studied (Pasho et al. 2011;Camarero et al. 2015Camarero et al. , 2018Kurz-Besson et al. 2016;Andivia et al. 2020;Vergarechea et al. 2021), showing that wood growth may be affected by such events. Conversely, other species are tolerant to dry conditions, such as holms and cork oaks, which are already being used to replace less tolerant species in high-elevation regions of NE Spain (Valladares et al. 2014). Cork oaks acquire protective mechanisms to buffer the typically large variability in the available water supply (Kurz-Besson et al. 2014), being nevertheless sensitivity to late spring precipitation.
Assessing the ecosystem vulnerability to drought events is a complex endeavour, which is influenced by the multitude of definitions and approaches to the problem. Therefore, and to tackle this issue, the Intergovernmental Panel on Climate Change (IPCC 2014) defined vulnerability as "the propensity or predisposition to be adversely affected, encompassing a variety of concepts including sensitivity or susceptibility to harm and lack of capacity to cope and adapt." Accordingly, drought vulnerability is driven by three components-exposure, sensitivity, and adaptive capacity (Murthy et al. 2015a;Alonso et al. 2019). The ecosystem's exposure to drought is assessed based in the extent, duration, and frequency of drought events in the region; sensitivity is related to the ecosystem's likelihood of being negatively impacted by drought; and adaptive capacity regards the ecosystem's ability to endure or recover (Murthy et al. 2015a;Engström et al. 2020). Vulnerability approaches have been used as a tool to mitigate adverse impacts on agriculture crops (Luers et al. 2003;Murthy et al. 2015b;Wu et al. 2017;Alonso et al. 2019) and forest ecosystems (Lindner et al. 2010;Sharma et al. 2015;Mildrexler et al. 2016;Swanston et al. 2018) in several areas of the globe. By developing a vulnerability approach applied to the Iberian Peninsula, forests in danger of suffering from wood losses may be highlighted, and mitigation measures contributing to minimize the impacts of these events may be locally promoted.
The objective of this work is to develop a vulnerability map of forests to drought over the Iberian Peninsula by taking advantage of drought indicators such as the Standardized Precipitation-Evapotranspiration Index (SPEI) and the remotely sensed Vegetation Health Index (VHI). This work aims at pointing into the direction of the most vulnerable regions, occupied by forests, to drought and warns for the need to implement forestation methods and policies that may mitigate species dieback. Here, we take advantage of a method already applied, with relevant results, for estimating maps of vulnerability focused on agricultural crops in India (Murthy et al. 2015b) and Portugal (Alonso et al. 2019) and adapt it to the forests of the Iberian Peninsula. The forest ecosystem is especially relevant here because of its strong socioeconomical impact in Portugal and Spain, and knowing where forests are more vulnerable to drought serves as a starting point to define adaptation strategies to climate change by policymakers and authorities. This study does not take into consideration tree species and specific biological traits, building on the derivation of a general vulnerability map developed from the climate perspective, in order to inform about where decision makers may direct their efforts. Further works may detail the biological characteristics of each tree species that may be found in these regions.

Study area and data
This work is focused on the Iberian Peninsula, located in the SW part of Europe. The IP is mainly characterized by three bioclimatic classes (Rivas-martínez et al. 2011;Gouveia et al. 2012Gouveia et al. , 2017: Mediterranean dry, Mediterranean oceanic, and temperate oceanic. Mediterranean dry climate is mainly located in the eastern IP, while the Mediterranean oceanic characterizes regions of western IP, most prominently over the SW corner of the peninsula. These regions are characterized by dry summers with propensity for intense, frequent, and durable droughts and mild, wet winters, which are responsible for the spring greenness peak . The N and NW IP are characterized by the temperate oceanic bioclimate, which has a distinct vegetation greenness annual cycle than its Mediterranean counterparts, with its peak in summer . The temperate oceanic is defined by a temperate, warm summer, and by no significant precipitation differences among seasons. Different bioclimatic regions show different associated ecosystems and types of forests (Barbati et al. 2007). The Mediterranean dry and oceanic regions are dominated by sclerophyllous vegetation, with mixed deciduous broadleaved forests interspersed in some regions, whereas the temperate oceanic region is characterized by deciduousand coniferous-broadleaved forests (Navascués et al. 2006;Barbati et al. 2007). Mediterranean forest species range from different pine trees, which are predominant in the region, to oaks and holms, and eucalyptus; furthermore, the temperate oceanic region is mainly composed of beech, pine, and eucalyptus forests (Navascués et al. 2006;Barbati et al. 2007;Uva et al. 2015;Nunes et al. 2020). The derivation of forest vulnerability to drought implies the use of several datasets in the development of its components (exposure, sensitivity, and adaptive capacity). The following subsections detail these datasets.

Standardized Precipitation-Evapotranspiration Index (SPEI)
The SPEI is a multi-scalar drought indicator that depends on precipitation and evapotranspiration information   (Harris et al. 2014), with a spatial resolution of 0.5 × 0.5°. SPEIbase was produced based on the FAO-56 Penman-Monteith estimation of potential evapotranspiration. The Penman-Monteith method is considered a superior method, being SPEIbase recommended for most uses including long-term climatological analysis . For this study, the 6and 12-month timescales in June and August are retrieved. The former represents early summer drought integrated with meteorological conditions in winter and spring, which includes the rainy season and the peak of vegetation greenness in Mediterranean bioclimates, while the latter represents an integrated picture of summer drought due to meteorological conditions within the hydrological year. SPEI was previously used by Alonso et al. (2019) to characterize exposure to droughts of agricultural crops. In contrast, Murthy et al. (2015a) described crop-generic exposure to drought by using seasonal and sowing period precipitation, which does not account for exposition to extreme temperatures nor changes in temperature patterns and variability due to climate change. Hence, by using SPEI both precipitation and temperature (in the form of evapotranspiration) ) are used to describe exposure.

Aridity index
Aridity is a natural phenomenon resulting from the imbalance between water availability and water demand, being a relevant element to characterize the climate of a given region with important implications for ecosystems, water management, and other human activities (Santos Pereira et al. 2009;Andrade and Corte-Real 2016). Aridity may be measured by the ratio between annual precipitation and evapotranspiration, which is commonly named as aridity index (AI). In this study, the AI was obtained from the CGIAR Consortium for Spatial Information (CGIAR-CSI) global aridity index and potential evapotranspiration climate database v2 (Zomer et al. 2006(Zomer et al. , 2008. This high-resolution (30 arc-seconds) AI was estimated by taking advantage of the Penman-Monteith equation to estimate reference evapotranspiration. The AI is classified as hyper arid ( AI < 0.03 ), arid ( 0.03 < AI < 0.20 ), semi-arid ( 0.20 < AI < 0.50 ), subhumid ( 0.50 < AI < 0.75 ), and humid ( AI > 0.75 ) (UNE-SCO 1979). The aridity index was also taken in consideration since it allows to classify a given region regarding its water balance under typical conditions of climate (Alonso et al. 2019).

Vegetation Health Index (VHI)
The VHI is a commonly used index to monitor drought based on remotely sensed in-formation (Kogan 1997;Bhuiyan et al. 2006;Kogan et al. 2012;Pei et al. 2018;Ribeiro et al. 2019). It consists of a linear combination between two components: the Vegetation Condition Index (VCI) and the Thermal Condition Index (TCI). The former takes advantage of information from the visible and near infrared portions of the electromagnetic spectrum, while the latter relies on the thermal infrared band. VCI is estimated based on the Normalized Differences Vegetation Index (NDVI) and aims at accounting for the vegetation water stress in each moment and pixel when compared to its maximum and minimum recorded values, and TCI is estimated with top-of-atmosphere brightness temperature (or land surface temperature, LST) and is used to assess the temperature stress of the vegetation also relating to the maximum and minimum recorded values. Recently, two works showed that the weights given to VCI (using NDVI) and TCI (using LST) may be equal or may depend on the limiting factor to vegetation growth (Bento et al. 2018(Bento et al. , 2020. In the present study, weekly values of VHI are retrieved from the NOAA center for satellite applications and research (STAR) at 4 km spatial resolution (https:// www. star. nesdis. noaa. gov). Here, this index is estimated as the simple average between VCI and TCI. Monthly data is currently available starting in 1982, and months from January to December are used in this study, to account for drought stressors to forests that may occur during the year.
With the objective of characterizing the sensitivity component, the VHI was used since it takes in consideration the vegetation state in a given month by means of the NDVI, which may be seen as an indicator of the susceptibility of the vegetation to drought (Kogan 1997). Other vegetation indices such as the Leaf Area Index could also be used here (Jump et al. 2017;Lecina-Diaz et al. 2021), but VHI presents the added value of also incorporating top-of-atmosphere brightness temperature, therefore being a more robust index to monitor drought in regions with more humid climates where vegetation growth is limited by temperature (Bento et al. 2018), and not only in water-limited vegetation that may be found in more semiarid regions (Bento et al. 2020).

Land cover
The European Space Agency (ESA) Climate Change Initiative (CCI) makes available global land cover maps at 300 m spatial resolution on an annual basis from 1992 to 2015 (http:// www. esa-landc over-cci. org/), describing a total of 22 land cover classes (2017). Furthermore, the Copernicus Climate Change Service (C3S) released land cover products for the years 2016, 2017, and 2018 (https:// cds. clima te. coper nicus. eu). These were developed to ensure continuity and thus to be consistent with the ESA-CCI land cover, being in the same resolution grid. The land cover used in this study merges the two datasets and spans the period between 1992 and 2018.

Adaptive capacity
As previously stated, this study aims at developing a vulnerability map based on climate features, without specifying tree species and biological traits. Therefore, the choice of indicators used to define adaptive capacity was made by taking in consideration previous works related to the physical and climatic characteristics of the region.

Elevation
Elevation in the Mediterranean forests, and thus in the Iberian Peninsula forests, plays an important role in the productivity of the forests, in particular concerning to temperature and moisture regimes associated with different insolation periods of the site topography. Elevation (among other features related with topography) affects the microclimates of the region, which in turn affects the distribution of vegetation and soil processes (Lozano-García et al. 2016;Helman et al. 2017;Molina-Venegas et al. 2017). With that said, elevation was included as an adaptative capacity component, since trees at different elevation levels will have different abilities to adapt to extreme climate events.
The elevation map used in this study comes from the National Oceanic and Atmospheric Administration's (NOAA) National Centers for Environmental Information (NCEI) (https:// ngdc. noaa. gov/ mgg/ topo/ globe. html). The dataset, named Global Land One-km Base Elevation Project (GLOBE), is a 30 arc-seconds longitude-latitude grid spaced Digital Elevation Model (DEM) (Hastings et al. 1999). Elevation is generally regarded as one of the most important adaptive capacity extrinsic indicators (Lindner et al. 2010;Lecina-Diaz et al. 2021).

Water table depth (WTD)
The relevance of groundwater in the adaptive capacity of vegetation in the Iberian Peninsula is widely recognized, with several studies mentioning its importance (Pinto et al. 2014;Kurz-Besson et al. 2016;Marques et al. 2019;Páscoa et al. 2020). A striking example is the study from Páscoa et al. (2020), where the authors developed a general method to identify regions with groundwater-dependent vegetation based on NDVI and validated the results by using a map of Water Table Depth (WTD). A global WTD (Fan et al. 2013) (https:// aquak now. jrc. ec. europa. eu) was used to infer the influence of shallow groundwater on forests. The WTD was estimated by taking advantage of several countries' government archives and available literature, whereas the missing data was filled by using a groundwater model forced by climate, terrain, and sea level developed by the authors (Fan et al. 2017). The map of WTD also has a spatial resolution of 30 arc-seconds. Furthermore, water table depth is seen here as a proxy to the rooting depth, i.e., water availability closer to the surface may help vegetation to efficiently recover after a drought event (Fan et al. 2017), thus being selected as an adaptive capacity indicator.

Solar irradiation
Another considered factor was solar irradiation. Nemani et al. (2003) described the geographic distribution of potential climatic constraints to vegetation growth derived from long-term climate statistics and concluded that vegetation growth is mainly limited by irradiation in the northern regions of the Iberian Peninsula. Solar radiation was also shown to be determinant for vegetation growth in other areas of the globe (Righi et al. 2007;Soares et al. 2009;Caron et al. 2018;Nardini et al. 2019).
Maps of global horizontal irradiation (kWh/m2) are obtained from Solargis (Solar resource map © 2019 Solargis; https:// solar gis. com) at a resolution of 30 arc-seconds. Global horizontal irradiation is the sum of direct and diffuse radiation received on a horizontal plane. This variable is used here since it is a limiting factor to vegetation growth in regions located in the N Iberia (Nemani et al. 2003).

Fire radiative energy
Forests located in the Mediterranean basin represent the most frequently fire-plagued regions (Pausas and Fernández-Muñoz 2012), with countries like Portugal (Trigo et al. 2006;Turco et al. 2019) and Spain (Rodrigues et al. 2019) being recurrently affected by these events, which in turn influence how vegetation recovers (Keeley 2009). Fire Radiative Energy (FRE) has been used to characterize fire intensity and therefore an indicator of extreme fire behavior that deeply affects forest resilience and difficult vegetation recovery after fires. Hence, FRE is chosen as an adaptive capacity indicator to map regions most vulnerable to wildfires.
Fire radiative power (FRP; megawatts) is obtained from the Satellite Application Facility for Land Surface Analysis (LSA SAF) (Trigo et al. 2011) FRP product. This pixelby-pixel variable consists of estimates of radiative power emitted by landscape fires from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) instrument on board the Meteosat Second Generation (MSG) series of geostationary satellites operated by the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT). The product is provided for the MSG disk every 15 min. A complete description of the product and its validation is available online at the LSA SAF site (http:// lsa-saf. eumet sat. int). The daily energy released by fire is estimated by integrating the FRP in one given pixel and one given day (Pinto et al. 2018). For data sampled every 15 min, the daily fire radiative energy (FRE) for each pixel p and day d is estimated as: where i represents the 15-min intervals within the day d , and the 0.9 factor converts the results into GJ. Daily FRE p,d is then accumulated for the period of study to obtain the map of FRE p . This variable is linked with fire intensity and fire severity, which are thus related to soil erosion and natural revegetation (Keeley 2009), being therefore used in this study as an adaptive capacity indicator, i.e., more FRE may handicap the ecosystem in its regeneration.

Data postprocessing
To pursue this study, all datasets must be geographically consistent and placed into overlapping grids. Thus, all datasets were spatially cropped to focus on the Iberian Peninsula, with longitudes from 10.75 W to 5.75 E, and latitudes from 35.75 N to 45.75 N. SPEI and VHI were retrieved for the common period 1982-2018, and land cover for the available period 1992-2018. Moreover, VHI and SPEI are used at the original weekly and monthly scales, respectively. Since all variables present spatial resolutions ≤ 4 km except the SPEI, the 4 × 4 km becomes the target resolution. As a result, the AI, elevation, WTD, FRE, and solar irradiation are interpolated to this grid using a nearest-neighbour method. To become comparable, the sparser 0.5° SPEI dataset is further downscaled into this resolution, by providing the same SPEI value to the analogous 4 × 4 km grid. Finally, the land cover map is converted into a forest cover map by giving each 4 × 4 km square the percentage of 300 × 300 m pixels associated with forest classes that are present within the square. Thus, according to ESA-CCI, the classes identified in this study as forest are tree covers with evergreen or deciduous broad or needle leaves, tree cover with mixed leaf type (broad and needle leaved), and mosaic tree and shrub with herbaceous cover. All the data processing was performed using Python 3.7 language.

Methods
The exposure of forests to drought is characterized by SPEI and the aridity index ( Table 1). The first indicates the exposure of a given region in a given month to drought considering different memory components (timescales) of the previous months' meteorological conditions, while the latter differentiates regions in terms of general water availability. The exposure component is then subdivided in eight indicators-seven related to SPEI and the remaining being the aridity index. Hence, the maximum and minimum 12-month SPEI in August are used since summer droughts are the most relevant to forest tree health (Pasho et al. 2011;Mildrexler et al. 2016;Schuldt et al. 2020). These indicators have information relative to the previous year conditions. The third indicator is the number of years in which the 12-month SPEI in August identifies extreme drought (SPEI ≤−2). The maximum and minimum 6-month SPEI in June are also considered. These indicators consider the occurrence of spring droughts, which are also a critical factor to the degradation of forests in the subsequent summer months (Kurz-Besson et al. 2016). The number of years in which the 6-month SPEI in June identifies moderate (− 1 ≥ SPEI ≥ −1.49) and extreme drought are the sixth and seventh indicators, respectively. SPEI thresholds for moderate and extreme droughts are according to Rhee and Cho (2016). The sensitivity component comprises five indicators, related to VHI and forest cover ( Table 1). The VHI metrics are considered annually and are divided in integrated and maximum indicators (Murthy et al. 2015a). The first indicator is the coefficient of variation (CV), i.e., the ratio of standard deviation to the mean, of accumulated VHI during the year, while the second is the CV of the maximum VHI of the year. Furthermore, to illustrate the frequency of droughts in 1 year, the third indicator is the CV of the number of weeks in the year with VHI lower than 20 indicating severe vegetation stress. The fourth indicator is the number of weeks in the total period of data in which VHI is less than 40, indicating vegetation stress. Finally, the fifth indicator of sensitivity is the mean forest cover for the available period .
Finally, the adaptive capacity component intends to describe the ability of a system to adapt to drought conditions (Lindner et al. 2010). This component should also be comprised of static parameters without intra-annual variation (Murthy et al. 2015b). Thus, the chosen indicators are the ones that may constrain the response of forest to dry conditions, namely elevation, water table depth, fire radiative energy, and annual horizontal irradiation (Table 1). Table 1 also shows the assumed functional relationships between indicators and the respective component, being either positive or negative. Lower values of maximum and minimum SPEI and AI represent more exposure to drought, and thus are assigned with a negative relation; more years with SPEI indicating moderate or severe drought represent increased exposure, being therefore assigned with a positive sign. Higher CV of accumulated and maximum VHI (lower mean VHI) and higher number of weeks with VHI lower than 40 indicate higher sensitivity; conversely, lower CV of the number of weeks with VHI less than 20 and lower mean forest cover (which may be related with transition zones from forested to not forested regions) indicate higher sensitivity, being the functional relationship negative. Finally, higher adaptive capacity is typically found on lower elevations, when water is available near the surface, when less energy from fires is released, and when less irradiation is present in summer, i.e., all variables present a negative functional relationship.
Exposure, sensitivity, and adaptive capacity components are obtained through principal component analysis (PCA) applied to the sets of variables described above (eight indicators to estimate exposure, five to sensitivity, and four to adaptive capacity). This method allows for a reduction of a dataset containing a large number of variables to a dataset containing fewer new variables, which are linear combinations of the original ones and represent the maximum possible amount of variability contained in the original data (a comprehensive description of the method may be found in Wilks (2011)). Furthermore, PCA allows to eliminate redundant information that may exist in the data due to correlations between variables (PCs are orthogonal among each other). Another advantage of using PCA is the ability to highlight spatial patterns within multivariate data (Abson et al. 2012). The indicators (Table 1) are further standardized and the PCA is applied to estimate exposure (based on eight indicators), sensitivity (based on five indicators), and adaptive capacity (based on four indicators). Once eigenvectors are found from the covariance matrix, the next step is to order them by principal components (PCs), highest to lowest, ordering the PCs in terms of significance. Here, the selection of the significant PCs was made by choosing the ones that represent more than 80% of explained variance. The inner product of the eigenvector matrix loadings ( L ) and scores ( S ) for those PCs is applied to obtain maps representing the vulnerability component ( v ): exposure, sensitivity, and adaptive capacity: where i is the pixel, j and N are the indicator and the total number of indicators of the vulnerability component v , respectively. Finally, n is the number of PCs used for a given vulnerability component, i.e., the number of PCs needed to achieve at least 80% of the explained variance. These are then transformed by multiplying by − 1 if the functional relationship obtained from the PCA is inverted, i.e., high exposure, high sensitivity, and high adaptive capacity are assigned positive signs. The procedure of PCA is then repeated, this time by using the transformed maps of exposure, sensitivity, and adaptive capacity as "indicators" to obtain forest vulnerability to drought (following Eq. 2). A simplified representation of the methodology is presented in Fig. 1.
Drought vulnerability values are then categorized in classes by dividing the obtained vulnerability map values into five intervals with the same probability of occurrence, in steps of 20%, i.e., dividing the empirical cumulative distribution function (ECDF) of the vulnerability map into five equal intervals. The five classes are then defined as less vulnerable (0-20% of probability of occurrence), moderately vulnerable (20-40%), vulnerable (40-60%), highly vulnerable (60-80%), and extremely vulnerable (80-100%) (Alonso et al. 2019).
Finally, forest types (according to ESA CCI land cover) that are found in the five vulnerability classes over the Iberian Peninsula are analyzed. Three different regions with large clusters of the extremely vulnerable class are selected and forest-types representative of these regions are detailed by verifying their frequency in those classes, i.e., by assessing the most vulnerable types of forest into these three regions.

Results
Principal component analysis is applied to indicators of exposure, sensitivity, and adaptive capacity. Figures 2, 3, and 4 show the explained variance of each PC for these three sets of indicators, respectively. The threshold established in this study to retain PCs is a cumulative explained variance of at least 80% of the original variance. The first five PCs are retained, accounting for ~ 86% of the variance for exposure (Fig. 2); for sensitivity ( Fig. 3) and for adaptive capacity (Fig. 4) the first three are selected, which explain circa 89% and 81% of the original variance, respectively.
The spatial patterns of eigenvector loadings for the retained PCs are displayed in Fig. 5. For exposure (Fig. 5,  top panel), the most relevant indicators for the PC1 (which These results show the importance of negative extremes of SPEI06 in June, being the most important factor to contribute into exposure from the first PC. Indeed, both the minimum SPEI06 in June (with negative sign), and the number of years SPEI06 in June indicate severe drought (with positive sign), show the larger contributions. The dichotomous sign is related to the fact that the variables are reversed, i.e., minimum SPEI and larger number of years in severe drought point to drought exposure. The importance of the first 6 months of the year in this region is paramount since it encapsulates the wet season and the peaks of greenness of the Mediterranean dry and Mediterranean and temperate oceanic bioclimates. It is also relevant to stress that the aridity index also has a large contribution, which is somehow expected, since the marked asymmetry between moist and dry regions in the Iberian Peninsula should imply different exposure characteristics when drought occurs. For sensitivity (Fig. 5, middle panel), PC1 (about 58% of the explained variance) shows positive contributions from the CV of VHI accumulated in the year, CV of maximum yearly VHI, and the number of weeks in the full period with VHI less than 40. The three indicators are similar in magnitude. Negative contributions to PC1 come from the CV of number of weeks in which VHI is less than 20 and the mean forest cover. Again, the signal of the loadings of the first PC is in agreement with the signal assigned to their functional relationship with sensitivity. PC2 (explaining about 21% of the explained variance) is characterized by positive contributions from all indicators except the number of weeks with VHI less than 40. PC3 (circa 11% of the explained variance) For adaptive capacity (Fig. 5, bottom panel), PC1 (which explains about 33% of the variance) has positive contributions from elevation, water table depth, and FRE (near zero), being the first two dominant indicators, and a negative contribution from horizontal irradiation. PC2 (25.3% of the explained variance) is characterized by positive contributions from FRE and negative contributions from the remaining indicators (with similar small magnitudes). PC3 (with an explained variance of 22.4%) presents a dominant positive contribution from horizontal irradiation, and similarly positive contributions from the remaining indicators. These results show very similar magnitudes between indicators in the first PC (except for FRE, which shows a small magnitude in the first PC, and is dominant in the second PC also with a positive sign), with elevation, FRE, and water table depth showing similar signs, which in turn were opposite to solar irradiation. This is on par with the above discussed, showing the predicament of forests to adapt when in high altitudes and when water is deeper. On the other hand, solar irradiation is an important asset for the development and recovery of forests (Nemani et al. 2003). However, the functional relationship of solar irradiation to adaptive capacity is set as having the same sign as the other indicators in Table 1, Fig. 2 Spatial distribution of the selected indicators to characterize exposure (numbered 1-8, see Table 1), and explained variance (%) of each principal component (PC). Percentages inside the graph indicate cumulative explained variance starting in the first PC which is not observed here. Indeed, solar irradiation has distinct roles regarding its importance at different times of the year, i.e., large amounts of irradiation may be beneficial to plant growth in winter, and prejudicial in summer. We used annual solar irradiation which does not allow to isolate the role played by this indicator, however, may be captured by PC3, which attributes equal signs to all variables.
The combination of the matrix eigenvector scores and loadings for the retained PCs results in the exposure, sensitivity, and adaptive capacity spatial patterns (Fig. 6). The map of exposure shows contrasting patterns between regions of the NW IP (which continues east along the coast until the Pyrenees), west IP (majority of the Portuguese territory), and SW IP, all showing negative values of exposure (less exposed) and regions of central towards the eastern coast of Spain, showing positive values (more exposed). Sensitivity shows the largest positive magnitudes (more sensitive) in the eastern Andalucía (W Murcia, S Castilla-la-Mancha), along the eastern coast of Spain, and in central Portugal. On the other hand, regions in the N and NW IP, the Portuguese coast from Lisbon towards N, and central Andalucía show the largest negative values of sensitivity (less sensitive). Finally, the adaptive capacity map displays negative values (less capable to adapt) in regions such as eastern Andalucía (W Murcia, S Castilla-la-Mancha), and central Portugal. Positive values (higher adaptability) are mainly located in the NW and N coasts of IP and parts of the eastern coast of Spain.
The principal component analysis is now applied to the newly determined exposure, sensitivity, and adaptive capacity maps, to obtain the vulnerability map (see Fig. 1). Since the threshold adopted in this study is 80%, the three PCs are used, accounting for 100% of the variance (Fig. 7a). Figure 7b shows that the first PC (explaining circa 48% of the variance) is characterized by a positive contribution from exposure and sensitivity, and negative from adaptive capacity, all showing similar magnitudes. This agrees with the functional relationships of vulnerability: more exposed,  Table 1) more sensitive, and with less capacity to adapt. Moreover, Fig. 7 c-e represents the spatial patterns of eigenvector scores, which shows that PC2 and PC3 contributions were not as relevant as PC1. Figure 8a displays the vulnerability map as obtained from the linear combination of the eigenvector scores and loadings of the three PCs.
The vulnerability classes are then defined by dividing the values of the vulnerability map in five equal-spaced intervals within its ECDF. The classes are then defined for cumulative probabilities as shown in Table 2, ranging from less vulnerable until extremely vulnerable. Hence, the map of vulnerability shown in Fig. 8a becomes the map of vulnerability classes shown in Fig. 8b, which further highlights regions with extremely vulnerable forests to drought.
The map shows clusters of extremely vulnerable forests to drought in areas such as the central region of Portugal, the eastern coast of Spain, and regions of the NE IP and central inland Spain. Three different regions, with large, forested areas, located in different regions of IP, and within different bioclimatic areas were highlighted in Fig. 8  The vulnerability map shows regions where the forest is distressed by drought events but does not specify the land cover which is most at risk. The four dominant forest types in the IP are the broad-leaved deciduous, needle-leaved evergreen, a blend of both, and the mix between mosaic trees and shrub. Figure 9 shows the percentage of these four types within each vulnerability class for the full area of study and the three specific highlighted regions. Results for the overall study region show a different behaviour between forest types, where needle-leaved evergreen and mosaic tree and shrub are found in larger number towards more vulnerable regions, broadleaf deciduous is mainly found in regions less vulnerable.
Forests situated in the central region of Portugal, in the limit between temperate and Mediterranean oceanic biozones (Region 1 in Figs. 8 and 9), show extreme vulnerability to droughts. Here, mosaic tree/shrub forests also represent a relevant land cover type in most of the vulnerability classes. However, needle-leaved evergreen forests are the ones with larger area within the extremely vulnerable class (~ 24%), closely followed by the mosaic tree/shrub (~ 20%).  Table 1) The largest percentage of mixed forests is also defined by extreme vulnerability (~ 3% of the area).
Region 2 (Figs. 8 and 9), located in the region of central-eastern IP, in the transition between the Mediterranean dry and oceanic bioclimates, shows another large cluster of vulnerability to drought in forested areas. The main forest types present in this region are mosaic tree and shrub, and needle-leaved evergreen land covers, with the first representing about 60% of the total forested area in the region, of which about 35% is extremely vulnerable, circa 15% are highly vulnerable, and 8% is vulnerable to drought events (accounting 58% of the 60% of this forest type in the three most vulnerable classes on the region; Fig. 9b). Needleleaved trees represent about 33% of the forests of the region, of which ~ 23% is in the extremely vulnerable class (Fig. 9b).
The dominant forest type in Region 3 (Figs. 8 and 9), located in the NE IP in Catalunya, is needle-leaved evergreen (about 72% of the total forested area).

Discussion
Forest ecosystems in the Iberian Peninsula are exposed to changes in climate that have already occurred and are expected to further intensify (Giorgi 2006;Naumann et al. 2018;Anderegg et al. 2020). In the current context of climate warming, drought events are more frequently observed, and present larger duration and increased intensity (Spinoni et al. 2015b(Spinoni et al. , 2018Páscoa et al. 2017Páscoa et al. , 2021. Forest productivity may be endangered since many tree species are drought sensitive (Galiano et al. 2010;Pessoa et al. 2014;Camarero et al. 2015Camarero et al. , 2018Gazol et al. 2017;Rubio-Cuadrado et al. 2018;Andivia et al. 2020), and even species tolerant to dry conditions may suffer from prolonged events. This opens the door to an ecological and economical problem in some regions of IP (Lopes 2013), with the migration of trees to other regions, forest dieback, wildfires, insect outbreaks, and a consequent change in the forest industry, which is an important economic driver for countries like Portugal and Spain (Rego et al. 2013;Pinho 2014;Uva et al. 2015;Sierra-Pérez et al. 2015). Vulnerability assessment and risk management approaches are important to tackle the effects of drought, highlighting the urgent need for metrics that can show where these changes are occurring. A growing number of studies have assessed the links between vegetation stress and mortality to drought and increased temperature (Jump et al. 2006;Kurz-Besson et al. 2014Sangüesa-Barreda et al. 2015;Anderson et al. 2018;Gazol et al. 2020), but the majority of them lack the inclusion of the three components of vulnerability (exposure, sensitivity, and adaptive capacity). In response, we propose the application of a robust statistical method which was previously applied to agricultural crops by Alonso et al. (2019). That study upgraded the method proposed by Murthy et al. 2015a, b, by automatizing the method by means of a PCA analysis of the vulnerability components. The use of PCA to determine forest vulnerability does not require the a priori knowledge of the functional relationships between the indicators to compute exposure, Fig. 5 Loadings of the chosen PCs of a exposure, b sensitivity, and c adaptive capacity. List of indicators: 1-maximum SPEI12 in Aug; 2-minimum SPEI12 in Aug; 3-number of years SPEI12 in Aug identifies severe drought; 4-maximum SPEI06 in Jun; 5-minimum SPEI06 in Jun; 6-number of years SPEI06 in Jun identifies moderate drought; 7-number of years SPEI06 in Jun identifies severe drought; 8-Aridity Index; 9-CV of annual sum of VHI; 10-CV of annual maximum VHI; 11-CV of number of weeks with VHI < 20; 12-number of weeks with VHI < 40; 13-mean forest cover; 14elevation; 15-water table depth; 16-fire radiative power; 17-horizontal irradiation sensitivity and adaptive capacity, which is required by the approach proposed by Murthy et al. 2015a, b. For this study, the PCA method was opted for its straightforwardness and to avoid the inclusion of errors in the relationships between the variables. Indeed, the present work is an adaptation of a previous study (Alonso et al. 2019) on cropping systems in Portugal. Similar variables are used to define exposure in both works since both SPEI and AI are variables only dependent  Following several studies focused on the recovery capacity of forests, we chose elevation (e.g., Lozano-García et al. 2016), solar irradiation (e.g., Nemani et al. 2003), and fire radiative energy (e.g., Turco et al. 2019) to complement WTD in this component. Other indicators were used in previous studies (Lindner et al. 2010;Lecina-Diaz et al. 2021), such as the regeneration characteristics, growth rate, or life expectancy of particular species, or forest management traits such as selective thinning (Vallejo et al. 2012). However, the goal of the present work was to develop a model that could be applied to the entire Iberian Peninsula focusing on the climate-driven impacts and not specifying tree species. The inclusion of tree-species and biological traits into the indicators would be a difficult endeavour to apply for such a large area and heterogeneous land occupation. Such analysis may be exploited in a future work focused on a smaller study area.
The application of PCA analysis to the three components of vulnerability shows that the most vulnerable forests to drought are those located in eastern coastal and interior Spain, and in central Portugal. These regions, although showing large vulnerability to drought, encompass different dominant species that should be taken into consideration.
In the case of Region 1 (central Portugal), the dominant needle-leaved species are pines (Pinus pinaster, Pinus pinea), eucalyptus, and the dominant broad-leaved deciduous is cork oak (Quercus suber) (Nunes et al. 2020). As previously mentioned, cork oak trees are characterized by a protective mechanism to buffer the potential unavailability of water, being nevertheless sensitive to late spring precipitation, which is decreasing and is expected to further decrease in the next decades (Pessoa et al. 2014;Kurz-Besson et al. 2014), resulting in the exacerbation of the occurrence of summer droughts that may result in the decline of these species. It is however important to refer that many of these stands are irrigated, since they embody a relevant factor to the economy of the country, being transformed in agroforestry sites (Aranda et al. 2007;Vessella et al. 2010;Camilo-Alves et al. 2020). Forests located in the inner central Portugal are commonly affected by wildfires, being one of the most fustigated regions of the IP in terms of burned  area (Trigo et al. 2016). Short-scale drought events are known to drive summer forest fire occurrences , which increase forest vulnerability to drought in that region. The extreme vulnerability of pine, eucalyptus, and cork oak forests to drought may have an important backlash in the Portuguese economy (Pessoa et al. 2014). Hardwood extracted from eucalyptus for grinding in the paper industry, softwood from pine trees (e.g., to be sawed), and cork extraction for the cork industry, represent an important source of national income (Lopes 2013;Rego et al. 2013;Pinho 2014;Uva et al. 2015).
In the case of Region 2, the combined effect of strong droughts with vulnerable tree species (Sánchez-Salguero et al. 2012), mountainous terrain, and water availability at deeper depths (Fig. 4), makes Region 2 woodlands extremely vulnerable to drought-induced forest dieback. Tree stands found in Region 2 are in its majority pine species (Sánchez-Salguero et al. 2012). Even Maritime pine (Pinus pinaster), a species adapted to dry conditions, showed vulnerability to droughts (Gouveia et al. 2012). Previous studies also showed the influence of drought on the radial growth of Black pine (Pinus nigra) (Martín-Benito et al. 2008;Candel-Pérez et al. 2012;Tíscar et al. 2018), Scots pine (Pinus sylvestris) (Candel-Pérez et al. 2012), and Aleppo pine (Pinus halepensis) (Gazol et al. 2017), concluding that these species are prone to decline due to droughts in this region. Furthermore, the number of drought events in central-eastern Spain has increased in the last decades and is expected to continue to increase in the future due to global warming (Spinoni et al. 2015a(Spinoni et al. , b, 2018Páscoa et al. 2017Páscoa et al. , 2021. Region 3, located in the Spanish autonomous region of Catalonia, is characterized by several different tree species, namely beech (Fagus sylvatica), holm oak (Quercus ilex), Scots pine (Pinus sylvestris), and Aleppo pine (Pinus halepensis) forests (Nunes et al. 2020). This area has experienced a gradual increase in temperature and evapotranspiration, and a decrease in precipitation in the last decades (Rubio-Cuadrado et al. 2018). Although located in a region with less horizontal irradiation (Fig. 4) and partly in the temperate oceanic bioclimatic region (Gouveia et al. 2012), the steep orography (Fig. 4) together with water stress due to prolonged drought events may disrupt local forest species (Lucas-Borja and Vacchiano 2018). The area of suitable habitat (for beech forests) is expected to decrease by 40 to 90% until 2070(Castaño-Santamaría et al. 2019. Another study (Chaparro et al. 2017) points to the vulnerability of beech forests in the same region, and to the ability of Aleppo pines to endure drought events. Furthermore, this study (Chaparro et al. 2017) also points to the sensitivity of broadleaved species (Holm oaks) to drought in Catalonia (Ogaya and Peñuelas 2006), which is also observed in Fig. 9. Indeed, Region 3 is the one with the steepest increase in the frequency of these species towards more vulnerable classes (although Region 2 also shows an increase, but with less area covered with this type of forest), contrasting with results from the overall IP. In the Pyrenees, drought-vulnerable species (e.g., Pinus sylvestris) have been declining due to both global warming induced severe and long droughts and by historical logging (Galiano et al. 2010;Camarero et al. 2011;Sangüesa-Barreda et al. 2015).
Apart from the three regions analyzed in detail, other regions present forests extremely vulnerable to droughts, such as the SW and NE Portugal, the SE coast of Spain, or the central-N Spain. The region of SW Portugal is recurrently affected by large wildfires, which threaten the health of the soil and destroy ecosystems and forests (Oliveira et al. 2021). The SE IP is mainly characterized by pinewoods (such as Pinus halepensis). These pine stands are expected to become increasingly more sensitive to intense dry spells, making them extremely vulnerable to drought (Gazol et al. 2017). On the other hand, less vulnerable regions are situated in central Andalucía, parts of central Spain, NW IP, or the Portuguese SE and Spanish Extremadura. Some of these are regions characterized by permanently irrigated lands used for agriculture, such as those located in the central Andalucía, Portugal SE or Extremadura, and the NW Spain near Zaragoza (https:// land. coper nicus. eu/ pan-europ ean/ corine-land-cover).
Few studies attempted to develop forest vulnerability maps for the regional scale (Neumann et al. 2017;Sommerfeld et al. 2018), without focusing on specific tree stands, and without assuming a priori knowledge to identify the functional relationships that link vulnerability and drivers. The inclusion of these a priori assumptions may amplify or dampen effects that may emerge at finer scales from interactions among multiple factors. Conversely, Forzieri et al. (2021) used a machine learning technique to assess the vulnerability of forests to fires, which IP results are somehow similar to those presented in this work, with the NW region less vulnerable, and regions such as central and SW Portugal, and NE IP being vulnerable to fires, and central-N Spain presenting the highest vulnerability to forest fires in the IP, which is in agreement with patterns observed in the vulnerability map developed in this work. However, it is important to stress that our work is not focused on vulnerability to fires, but instead on the vulnerability of forests to drought occurrences (still, forests in regions severely burnt will need more time to recover). Furthermore, the methodology described in this study may be applied to forests located in other parts of the world, always being careful with the months and timescales chosen to capture exposure to SPEI. Both SPEI and VHI are globally available, making them perfect indicators to use on other parts of the world. Regarding adaptive capacity, many other indicators may be chosen, depending on the region of interest. In a plane study area focused on agroforestry, one can replace elevation by (e.g.) insect outbreaks. Moreover, this work uses information, encompassing the years 1982-2018, and intends to design a vulnerability map that is contemporary to this period. However, this study opens the door to future works focused on climate change by using the output from ensembles of regional climate models, evolving vegetation maps, and probabilities of ignition based on fire risk maps, that could be used to estimate forest vulnerability in different scenarios of greenhouse gas emissions.
The forest vulnerability map developed in this study highlights the general dependence of forest growth and therefore production of forest goods, to drought episodes. However, it is worth stressing that the vulnerability assessment here presented is based on droughts that occurred in the last decades. It appears inevitable that the Iberian Peninsula forests will have to cope with significant extreme drought events in the near future (Spinoni et al. 2018), making this a pressing matter that needs to be addressed by decision-makers and authorities responsible for the conservation of these forests. Indeed, forests play a fundamental role in the socioeconomic welfare and in the sustainable development of countries like Portugal and Spain. It is also relevant to emphasize that this study does not account for changes in forest cover in the future and should be viewed as a tool to understand forest vulnerability. Also, the vulnerability assessment performed here only takes in consideration responses to drought episodes, which are important but not the only risk factor to forests and consequently timber harvesting. Indeed, other mechanisms such as wildfires (Westerling et al. 2006;Schoennagel et al. 2017), insect infestations (Roy et al. 2014;Ramsfield et al. 2016;Potter et al. 2019), or changes in land use (Stork et al. 2009;Rathore et al. 2019) also contribute to degrade the capacity of forests to endure through episodes of drought. Forthcoming works may take advantage of the methodology to further analyze forest vulnerability into the future. Droughts, together with wildfires, insect outbreaks, and changes in the land use have already destroyed several forests worldwide, causing ecological, economical, and social costs, and further releasing large quantities of carbon dioxide, increasing feedbacks towards more global warming.

Conclusions
The sustainable exploration of forests is a key economic driver for countries like Portugal and Spain, located in the Iberian Peninsula, in the westernmost tip of the Mediterranean region. Thus, the adequate management of forests in the IP must be seen as a focal point by governments and specialized authorities. Such management takes new proportions due to the recurrent episodes of forest decline reported in the region, which may be explained by anthropogenic air pollution, pathogens, and climatic factors. Concerning the latter, drought episodes are typically identified as the most common explanation for the massive dieback of forests observed in the IP. Hence, it is paramount to assess forest vulnerability in the region.
In this work, and following the definition of vulnerability from the IPCC, the exposure, sensitivity, and adaptive capacity of forests to drought are addressed. By combining them, a vulnerability map is obtained. Results show that: • The most vulnerable forests to drought are those located in central Portugal, in the central-E Spain, and in the NE peninsula. Other examples of regions with vulnerable forests to drought are those located in the SW and NE Portugal, and in the region of central IP. • Relevant indicators to exposure are those pointing to drought in early summer accounting for the conditions of precipitation and evapotranspiration of the first semester of the year, encapsulating the rainy season and the peaks of greenness characteristic of the forests of IP. The annual balance between precipitation and evapotranspiration, measured as the aridity index, is also an important indicator of exposure to drought. Furthermore, VHI and the percentage of forest cover are pertinent indicators to assess the sensitivity of forests to drought. Finally, elevation, water table depth, fire radiative energy, and solar irradiation are important indicators to assess the adaptive capacity. • Mosaic tree and shrub is the forested land cover that is more affected by extreme vulnerability to droughts in the Iberian Peninsula, followed by needle-leaved forests, dominant in central Portugal, and NE IP.
These results may be used to apply forest management policies in regions endangered by forest decline due to recurrent drought episodes. Typical practices used to manage forests, such as thinning, cleaning, reforestation, and other forestry techniques may need to be adapted to further account for changes that may arise due to prolonged, more frequent, and more severe drought. Henceforth, wildfire and fuel management may also need to be revisited and further adapted.
Author contribution Data extraction, development of statistical models, and writing-original draft preparation, V.A.B.; Conceptualization, A.R. and C.M.G; writing-review, V.A.B., A.R., I.V., and C.M.G. All authors have read and agreed to the published version of the manuscript.
Code availability Not applicable.

Declarations
Ethics approval Not applicable.

Consent to participate
All authors consent to participate into the study.

Consent for publication
All authors consent to publish the study in a journal article.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.