Higher site productivity and stand age enhance forest susceptibility to drought-induced mortality

Warmer and drier conditions increase forest mortality worldwide. At the same time, nitrogen deposition, longer growing seasons and higher atmospheric CO 2 concentrations may increase site productivity accelerating forest growth. However, tree physiological studies suggest that increased site productivity can also have adverse effects, reducing adaptation to drought. Understanding such intricate interactions that might foster tree mortality is essential for designing activities and policies aimed at preserving forests and the ecosystem services they provide. This study shows how site factors and stand features affect the susceptibility of Scots pine to drought-induced stand-level mortality. We use extensive forest data covering 750,000 ha, including 47,450 managed Scots pine stands, of which 2,547 were affected by mortality during the drought in 2015 – 2019. We found that the oldest and most dense stands growing on the most productive sites showed the highest susceptibility to enhanced mortality during drought. Our findings suggest that increasing site productivity may accelerate the intensity and prevalence of drought-induced forest mortality. Therefore, climate change may increase mortality, particularly in old and high-productive forests. Such exacerbated susceptibility to mortality should be considered in forest carbon sink projections, forest management, and policies designed to increase resilience and protect forest ecosystems.


Introduction
Warmer and drier conditions associated with climate change result in accelerating forest mortality worldwide (Allen et al., 2010;Millar and Stephenson, 2015;Seidl et al., 2017). In Europe, the frequency and intensity of droughts have increased during the last 30 years (Spinoni et al., 2018). Especially the recent sequential droughts in 2015-2019 initiated large-scale tree mortality in central European forests (Brun et al., 2020;Schuldt et al., 2020). With ongoing climate change and an accumulation of stress conditions, many regions are expected to experience increasing drought-induced forest mortality (Sánchez-Pinillos et al., 2021).
But not only increasing droughts significantly affect forests. Higher temperatures, increasing atmospheric nitrogen deposition, and increasing CO 2 concentration alter forest ecosystems functioning, mostly increasing site productivity and stimulating forest growth (Lindner et al., 2014;Pretzsch et al., 2014b;Reyer et al., 2017). A commonly used indirect proxy of forest site productivity is the site index (SI) estimated using stand height at a specific age, which reflects the overall effect of all environmental factors on tree growth (Bravo-Oviedo et al., 2010). SI is used as a measure of growth conditions (Pau et al., 2022;Socha et al., 2021a). In temperate European forests, increased site productivity and accelerated growth have been documented since the last century (Appiah Mensah et al., 2021;Bontemps and Bouriaud, 2014;Kahle et al., 2008;Pretzsch et al., 2014b;Socha et al., 2021a). This increased growth has translated into increasing sequestration of anthropogenic carbon emissions (Le Quéré et al., 2018;Nabuurs et al., 2013). Yet, there are also first signs of carbon sink saturation (Nabuurs et al., 2013;Xu et al., 2021), and accelerated growth can cause an imbalance between growth and mortality rates resulting in a decrease in forest biomass stocks (Brienen et al., 2020). Increased site productivity and faster growth shorten tree lifespan and increase the rate of forest mortality (Brienen et al., 2020). Moreover, accelerated growth supports trees to reach a certain height in a shorter time and larger trees are most prone to drought-induced mortality (Bennett et al., 2015;Bigler, 2016;Merlin et al., 2015;Stovall et al., 2019). This phenomenon is linked to tree physiology. Fast-growing, taller plants have wider water-conducting conduits, which during drought are more vulnerable to xylem embolisms, hydraulic failure and carbohydrate starvation (Adams et al., 2017;Olson et al., 2018;Choat et al., 2018). However, there are also studies suggesting that the relationship between tree size and intraspecies mortality may be negative (Trugman et al., 2021, Ripullone et al., 2020. This may be because taller trees have a number of adaptations, including deeper root systems that allow better access to water during drought, more efficient water use and transport, and greater water uptake and storage capacity, which mitigate the decrease in water potential associated with path length (Fernández-de-Uña et al., 2023). Thus, the links between tree height and drought-induced mortality are still unclear and the underlying physiological mechanisms are not well understood and require further research particularly at the intraspecific level (Fernández-de-Uña et al., 2023).
Trees can adapt to dry site conditions (Zanetti et al., 2015), this capacity is particularly important in the case of long-term changes in soil moisture and the negative trend in soil moisture anomalies observed over recent decades in Europe (EEA, 2023). However, adaptation is connected, inter alia, to the plasticity of the root systems, which may decrease with age (Correa et al., 2019), which in turn may result in greater drought vulnerability of larger and also older trees. Moreover, in older and higher trees, photoassimilates and water must be transported over longer distances so that the balance of non-structural carbohydrates along the tree may be disturbed during drought (Adams et al., 2017;Hesse et al., 2021;Sevanto, 2018). Fine roots of large mature trees are particularly affected by non-structural carbohydrates deficiency during drought because they are the furthest away from the leaves (Landhäusser and Lieffers, 2012;Ryan et al., 2014). Therefore, drought resilience may decrease with tree age and height (Martínez-Vilalta et al., 2012).
Consequently, higher and older stands growing on more productive sites may show increased susceptibility to drought and increased canopy mortality. However, evidence demonstrating increased drought-induced mortality of the oldest, most productive stands is still lacking (Brienen et al., 2020;Choat et al., 2018). Moreover, many studies on tree mortality mainly consider the impact of drought and other exogenous factors Choat et al., 2018), but spatial patterns of mortality and underlying factors remain difficult to understand (Changenet et al., 2021). Therefore, from the point of view of the possibility of preventing exacerbated forest mortality, it is crucial to identify the much less well recognised endogenous factors such as stand and tree characteristics (McNellis et al., 2021). Stand height, density, volume or stand age can be useful in determining vulnerability to mortality. However, most studies related to drought-induced mortality are conducted at the individual tree or plot level (Allen et al., 2010;Allen et al., 2015;Choat et al., 2018). Nevertheless, the mortality process is characterised by high spatial variability, making it difficult to identify the determinants of this phenomenon over large areas based on the individual tree and plot data (Allen et al., 2010). For example, research on aspen mortality in western Canada has shown that patchiness hinders a reliable assessment of the mortality process and biomass losses across the region solely by extrapolating from a network of monitoring plots (Hogg et al., 2008;Michaelian et al., 2011). To overcome the problems associated with the high spatial variability of the mortality process, it seems most reasonable to combine local-and landscape-scale estimates of forest mortality and their underlying causes, as suggested by Hartmann et al. (2018). This highlights the need to develop new approaches to model mortality patterns at different spatial and temporal scales. Therefore, the use of wall-to-wall forest mortality data from across the study area can improve understanding of the processes underlying forest mortality, which is crucial for predicting the effects of changing climate conditions on forest ecosystem dynamics, in particular CO 2 sequestration and other ecosystem services (Bottero et al., 2021;IPCC, 2018;Perring et al., 2015).
To address existing challenges, here we predict the stand-level risk of drought-induced mortality and explain how stand features, site factors and site productivity affect the susceptibility of Scots pine stands to drought-induced mortality in Central Europe. Instead of individual tree or the sample-plots observations most frequently used to study mortality patterns, we used wall-to-wall data on forest mortality covering 750,000 ha, including 47,450 managed Scots pine stands from southern Poland, of which 2547 were heavily affected by mortality during the drought in 2015-2019. By drought-induced mortality we also mean mortality by other causes facilitated by drought such as insect and fungal pathogen gradations (Boczoń et al., 2021;Mezei et al., 2017;Nordkvist et al., 2023). To eliminate the impact of background mortality and therefore avoid interpreting self-thinning and random events involving small gaps in stands as drought-induced mortality, we defined stand-level mortality as those in which the area occupied by dead trees removed within sanitation felling exceeded assumed threshold of 1000 m 2 . We have focused our research on Scots pine, which is the most widely distributed tree species on Earth (Eckenwalder, 2009) and one of the ecologically and economically most important tree species in Europe (Buras et al., 2018;Buras and Menzel, 2019). In Poland, Scots pine is the main forest forming tree species occupying over 60% of the forest area (Pasichnyk et al., 2022).
The aim of the first stage of the analyses was to determine how the intensity of mortality in individual years expressed by volume of dead trees is affected by drought. The aim of the second stage of analysis, was to indicate stand and site characteristics determining the risk of standlevel mortality during the period of sequential droughts in 2015-2019. We hypothesised that (h1) intensity of Scots pine mortality is driven by climatic water balance (CWB), moreover we hypothesised that the risk of mortality during drought is increased in the oldest Scots pine stands (h2) growing on most productive sites (h3). We also assumed that susceptibility to drought might be affected by other forest and site features (h4), particularly stand density, which may exacerbate the adverse effects of drought by increasing competition between individual trees.

Methods
The study area is located in southern Poland. Of the total area of 2.5 milions ha, forests occupy approx. 30%, i.e. ~750,000 ha ( Fig. 1). Scots pine (Pinus sylvestris), with a share of 66% is the dominant tree species in the studied region (Fig. 1b). The Scots pine-dominated stands (with the share of Scots pine more than 70%) included in the study are subject to planned silviculture and are managed by State Forests.

Sanitation felling data
Sanitation felling data were used to analyse the effect of drought on mortality intensity of Scots pine in the research area and for modelling stand-level mortality probability. Sanitation felling is a harvest of trees for the purpose of removing insects or diseases from a stand (Pietzsch et al., 2021). In forests managed by State Forests in Poland sanitation felling are mandatory for forest protection reasons. Depending on the occurrence of tree mortality, sanitation felling is carried out on an ongoing basis. The analyses used data on the total volume of timber harvested within sanitation felling in each stand in each year over the last decade. Thanks to the annual resolution, the sanitation felling data allows to follow the pattern of forest mortality at stand-level scale in different years.
For stands with mortality in which sanitation felling is carried out, evidence cards are established in State Forest information system. The cause of mortality and volume of the dead trees are documented. The volume of each dead tree removed as sanitation felling is individually measured, however the sanitation felling database contains the standlevel records with aggregated volume of dead trees removed each year attributed to one of two groups of causes of mortality: 1) drought/pests and 2) windbreaks. In this study only the data from group 1) -drought/ pests, were used.
In the majority of the stands subjected to sanitation felling, the total volume of trees harvested in sanitation felling did not exceed 2 m 3 ( Supplementary Fig. 1). In all age classes, the size of total volume of trees removed from single stand within sanitation felling in 2015-2019 most often did not exceed 2% of the stand volume ( Supplementary Fig. 2). Mortality was present in all age classes, but the number of stands in which it was found was lowest in the 20-40 year class, in which sanitation felling were undertaken in 5.3% of stands ( Supplementary Fig. 2). In subsequent age classes, the share of stands with sanitation felling was gradually increasing from 12.3% in the 40-60 age class to 37.2% in stands aged 120-140 years ( Supplementary Fig. 2).
In intensively managed Scots pine stands growing in the research area, the phenomenon of tree self-thinning is greatly reduced by silvicultural treatments repeated in each stand every 5-10 years, which reduce competition and considerably reduces the background mortality. Despite this, individual cases of self-thinning are observed in managed stands. Therefore to further eliminate the impact of background mortality and therefore avoid interpreting self-thinning and random events involving small gaps in stands as mortality caused by drought and pests, we defined stands affected by mortality as those in which the area of sanitation felling exceeded 1000 m 2 . The predefined threshold does not refer to the amount of dead wood removed, which prevents overestimating mortality probability for older and denser stands with higher SI.

Relationship between intensity of mortality and climatic water balance
To determine the relationship between tree mortality intensity and drought, we analysed the relationship between the yearly aggregated volume of dead trees harvested on the research area as sanitation felling and climatic water balance (CWB) in particular years. Short-term (bimonthly) data on climatic water balance (CWB), which is calculated as a difference between precipitation and potential evapotranspiration are publicly available as part of the Agricultural Drought Monitoring System program, developed by the Institute of Soil Science and Plant Cultivation -National Research Institute in Puławy (MINROL, 2023). Bimonthly CWB is published with a 10-day step. To assess the intensity of drought in individual years, averaged CWB values were determined for the period April-September and for the period May-August. We chose these two overlapping periods to determine whether CWB related to the growing season or CWB related to the shorter season with higher temperatures was more relevant for Scots pine mortality. To consider the cumulative effect of drought in subsequent years, apart from the direct correlation of mortality intensity expressed by the volume of dead trees with CWB in a given year, we also analysed the lagged effect of drought. For this purpose, we correlated the mortality in a given year with the averaged CWB for the period April-September and May-August of the previous year. A linear regression model was used in this analysis. To check autocorrelation of residuals we used Durbin-Watson test. To meet our objective the mortality data from the period 2009-2019, attributed as drought/pests were aggregated into individual calendar years for the entire study area.

Development of the stand-level mortality probability model
The aim of the second stage of analysis was to develop a stand level mortality probability model that describes the risk of mortality in a given stand. In order to identify stand and site factors that predisposed stands to mortality during drought, we performed analyses for the severe drought period 2015-2019 and the period 2009-2014, in which the drought occurred only accidentally in low intensity. The analyses were carried out according to the scheme presented in Fig. 2, which include: 1. Data acquisition. 2. Determination of the variables used for modelling. 3. Model development and testing. 4. Prediction of the mortality probability.

Fig. 1.
Location of the study area within the Scots pine distribution range in Europe (a). Classification of forests in the study area into Scots pine-dominated stands (the share of Scots pine more than 70%) and other stands (b).

Measurement of the detailed characteristics of individual stands using airborne laser scanning point clouds
To measure characteristics of individual forest stands, we used airborne laser scanning (ALS) data. The ALS point clouds were collected between 2011 and 2015. The density of the point cloud was at least four pulses per square metre. The data were collected in leaf-off conditions (October to April) with the maximum scan angle at a level of 25 • . The ALS data were processed using algorithms available in lidR package for R (Roussel et al., 2020): the point clouds were normalised to the height above ground using the "normalize_height" algorithm; for creating the Canopy Height Model (CHM) we used the "p2r" algorithm which uses the maximum point height as value for CHM (the subcircle parameter was set to 0.3 m); finally the single tree detection was performed using the "lmf" algorithm based on a local maximum filter with varying search window. For ALS data analyses, the area occupied by Scots pine stands was divided into 20 × 20 m grid cells. Based on ALS point clouds, tree density (trees/ha) and top height (m) for each 20 × 20 m grid cell were calculated to characterise stand parameters. The top height for each 20 × 20 m grid cell was calculated as the mean value of the maximum values of CHM obtained from the four 10 × 10 sub-cells, according to the approach proposed by Socha et al. (2020b). According to a previous study conducted for Scots pine stands in Poland, the automatic detection of trees based on ALS data works well when trees higher or equal to 60% of the top height of the stand are considered (Socha et al., 2020a). Thus, for calculating the tree density (trees/ha), we used only the detected treetops with a height equal to or higher than 60% of the top height.

Determination of remaining stand and site characteristics
ALS-derived top height and stand age retrieved from the Forest Data Bank (FDB, 2018) were used to calculate the SI, which was used to describe site productivity. The study concerns managed Scots pine stands that are even-aged and mainly originate from artificial regeneration. As management plans have been in place for nearly 150 years in the study area, the age of the managed Scots pine stands is well documented. SI was calculated using the base age 100 years according to the model developed for Scots pine in Poland (Socha et al., 2020c). The remaining stand and site properties were estimated using the data from FDB: age, dominant species, percentage share of dominant species, soil type, site type, the humidity of site (wet, fresh, dry). The humidity of site was assessed by experts during habitat mapping, based on analysis of the soil, topography and species composition of undergrowth plants. In addition, the modelling used the average May-August CWB values for the drought period 2015-2019 and the period 2009-2014 assigned to each stand and tested the suitability of elevation determined from ALS data.

Verification of mortality events in individual stands
The Scots pine stands indicated in the database as affected by mortality in 2015-2019 attributed to drought/pests have been verified using the land cover change masque obtained from the Landsat imagery. The first step of this process consisted of two land cover classifications using imagery from July 29, 2013 (Landsat 5) and August 31, 2019 (Landsat 8). Land cover classes included deciduous forest, coniferous forest, low vegetation, and other areas. Reference samples for training and validation were manually delineated based on the Landsat images; 275 samples were prepared for 2013 and 293 for 2019. The classification was performed using the Random Forest (RF) algorithm. RF is a nonparametric, ensemble classifier consisting of a combination of decision trees (Breiman L., 2001). Reference samples for classification were randomly divided into training and validation samples in a 70:30 ratio. The obtained land cover maps achieved an overall accuracy of 97.4% and 95.8%, respectively, for 2013 and 2019. The producer's and user's accuracy for the coniferous forests class, which was further analysed, were 100% and 97.3% for 2013, and 97.1% and 96.1% for 2019. The coniferous forests masque derived from land cover classifications was overlaid to acquire coniferous forests change masque in the next step. The change masque was then further visually inspected for errors, such as errors related to residual cloudiness/shading. The coniferous forest change was then spatially intersected with the stand polygons from the FDB, which enabled us to verify the occurrence of mortality in each of the Scots pine stands indicated in the sanitation felling database. Stands with sanitation felling caused by drought or pests were assigned as class 1 and stands without fellings were assigned as class 0. In the first variant of the model class 1 was limited to stands in which the area of forest cover loss due to sanitation fillings was larger than 0.1 ha. The adoption of the solution with a minimum area instead of a minimum volume was justified by the fact that if the same proportion of trees dies in terms of volume or number, the probability of mortality in all stands, regardless of their volume and density, will be the same. In order to test whether the introduction of a minimum sanitation felling area (0.1 ha) has an effect on the results, we also developed a model in which all stands included in the sanitation felling database were included without minimum area threshold. In total, 2547 Scots pine stands with sanitation felling in 2015-2019 assigned to class 1, and 44,903 stands assigned to class 0 were used to modelling stand-level mortality probability.

Model fitting
The model of stand-level mortality probability was developed using the logistic model from the family of the Generalised Additive Models (GAMs; Hastie and Tibshirani, 1990). As our aim was to estimate the risk of occurrence of mortality in the stand, we decided to use logistic regression. This was justified on the grounds that direct modelling of the proportion of dead trees at the individual stand level can be problematic for several reasons. Firstly, the approach taking into account all deaths of individual trees could be distorted by background mortality. Self-thinning and random events resulting in small gaps in stands could then be interpreted as mortality caused by drought and pests. Using a threshold in the form of a minimum volume of dead trees, on the other hand, could result in the mortality risk being overestimated in stands with the highest volume, i.e. older stands and those growing on more fertile sites. Based on the literature review (Woolley et al., 2012;Alenius et al., 2003;Salas-Eljatib et al., 2018), we concluded that a logistic model approach would be most appropriate for our study. However, the referred studies were mostly focused on modelling the mortality of individual trees. The GAM is a semi-parametric extension of the generalised linear models (Hastie and Tibshirani, 1990) and allows establishing a direct relationship between the explained variable with the individual explanatory variables holding all other variables fixed (Larsen, 2015). The general form of the GAM model is described by the following equation.
where E(Y) is the expected value of the response variable Y with an appropriate distribution, f j is a smooth function of the covariate x j , β 0 is an intercept term, and g is the link function (Pedersen et al., 2019). In case of this study the binomial distribution of response variable and logit link function were used. Two classes (0 and 1) were defined for logistic regression modelling: class 1 -stands with mortality, and class 0 -stands without mortality (Fig. 3). The stand-level mortality probability was calculated for the entire drought period 2015-2019 and nondrought period 2009-2014.
During model development, we tested the influence of site (CWB, elevation, soil type, site type, aspect, the humidity of site -wet, fresh, dry) and stand (SI, top height, stand density, age) characteristics on mortality probability. Redundancy of independent variables was estimated using correlation with other independent variables (Dormann et al., 2013). Finally, we selected only statistically significant and not redundant site and stand variables that affect the mortality probability. The model development was performed using the R package "mgcv".
Smoothness parameters were chosen using restricted maximum likelihood selection (Wood, 2011). The final model contains five most important independent variables: stand age, SI, stand density, aspect, and CWB. We used thin plate regression spline ("tp") as smooth function for four predictors (age, SI, stand density, CWB) and cyclic cubic regression splines ("cc") for terrain aspect. Ten basis functions were used ("k = 10″) for stand age, SI and aspect, and five basis functions ("k = 5″) for stand density and CWB.
To interpret the influence of the analysed predictor variables on the probability of stand-level mortality we used the partial effect plots. The partial effect plots were created separately for each predictor (stand age, SI, stand density, aspect, CWB) by performing model predictions with keeping values of other variables at their average value. The effect of the SI on the probability of mortality during a drought is most significant in older stands. Therefore, an additional partial effect plot was created to better illustrate the effect of SI on mortality probability for stands in age 50 and 100 years (+/− 10 years) and average CWB in May-August period equal − 20 and − 130 mm. To do that, the mean density and aspect were calculated for stands of age from 40 to 60 years and from 90 to 110 years, respectively, and model predictions were made for the whole range of SI observed in the input data.
For creating a map of stand-level mortality probability we developed an additional GAM model where we incorporated the spatial coordinates X and Y as predictor variables (coordinate system: EPSG 2180). Main effects with interaction of both coordinates were used in the model. It was assumed that accounting for the spatial structure of the data could help to explain more variability of mortality and provide better results when predicting the probability of mortality for the whole study area. Based on that spatial model the final map of stand-level mortality probability for the Scots pine stand was created.
There is no consensus in the scientific literature on the optimal approach for developing predictive models using the logistic regression approach in cases of severe class imbalance (Crone and Finlay, 2012;Salas-Eljatib et al., 2018). We, therefore, compared two approaches: i) fitting the model with all available observations from 47,450 stands with the class imbalance (44,903 standsclass 0, 2547 -class 1) and ii) fitting the model with random subsampling of class 0 (no mortality) to 2547 observations so that the training data are balanced (50%/50%, Supplementary Fig. 3). The analysis was performed for models including spatial component (X, Y coordinates). As there were no significant differences and slightly better validation statistics were observed in the case of the unbalanced dataset (Supplementary Table 1), the final version of the mortality model was developed on imbalanced data. Model validation was performed using the bootstrap approach with 500 iterations. In each bootstrap iteration, a sample with replacement from the input dataset was taken for model training. Then, prediction for the hold-out set was performed, and the performance metrics were calculated based on hold-out observations. On average, 37% of observations were used for validation at each iteration.

Relationship between mortality and climatic water balance
The intensity of Scots pine mortality is driven by CWB, however, the period for which the CWB is calculated plays a crucial role. The intensity of mortality in a given year expressed by the total volume of removed dead trees was not significantly affected by the averaged CWB in the period April-September. The CWB during May-August supposed to represent growing season conditions, explained about 46% of mortality intensity (p = 0.0215). In contrast, we found a much stronger lagged effect of drought on the mortality. The intensity of mortality was stronger correlated with CWB in the April-September period of the previous year (R 2 = 0.7470, p < 0.0001). However, the intensity of mortality was stronger driven by CWB in the May-August period of the previous year, which explains 89.5% of its variation (R 2 = 0.8951, p < 0.0001; Fig. 4).
The results showed that in Scots pine stands from the research area drought intensification occurred in 2015-2019, when CWB in the period May-August decreased under about − 120 mm, while lagged increased mortality was observed in 2016-2019 (Fig. 5).

Effect of site and stand factors on mortality probability during drought
With the exception of elevation, which shows strong concurvity with CWB, all other analysed independent variables were found to be significant predictors of probability of stand-level mortality during the drought period 2015-2019 (Table 1) and were used in the model. The explained deviance of the model was 31.0% (n = 47,450). In contrast to the model for the drought period, the model developed for the nondrought period (2009-2014) has a very low explained deviance (5.3%). This means that, apart from the drought period, the stand and site characteristics analysed are not significant predictors of mortality probability.
We showed that stand age, site productivity expressed by SI, and stand density are significant factors predisposing Scots pine stands to mortality within the drought period 2015-2019 (Fig. 6). Stand age affects the probability of mortality to the greatest extent (Fig. 6). The probability of mortality was lowest for stands <75 years, Fig. 6), and from the age of 75 to around 120 years the probability of mortality rapidly increased. The next factor significantly affecting mortality probability was the site productivity expressed by the SI. The mortality probability was lowest in stands with SI below 20 m (Fig. 6) and systematically increases with increasing SI. However, the effect of SI is agedependent and is particularly significant in older stands (Fig. 7d, Fig. 8). Stand density was also an significant determinant of mortality probability. An increase in stand density to about 450 trees per ha resulted in a systematic increase in mortality probability (Fig. 6). Probability of mortality was significantly correlated with CWB. Similar to the 2009-2019 time series, the analyses for individual stands for the period of 2015-2019 also found that mortality significantly increased when CWB in May-August was below − 120 mm (Fig. 6).
Increased probability of mortality was also found on the warmest and the driest south-western exposures (aspect 202.5-247.5 • , Fig. 6). On north-facing slopes (aspect 337.5-22.5 • ), the mortality probability was  The effect of the SI on mortality depends on the age of the stand and the average CWB in May-August. For 50-year-old stands with a CWB of − 20 (Fig. 7a), the modelled mortality probabilities were close to 0. For 50-year-old stands with a CWB of − 130 (Fig. 7b) and 100-year-old stands with a CWB of − 20 (Fig. 7c), there was a slight risk of mortality in stands growing on the most fertile sites. A very strong effect of SI on mortality probability was observed for 100-year-old stands with average CWB of − 130 (Fig. 7d). It can be observed that for oldest stands (age = 100) with strongly negative CWB (− 130 mm) with SI above approximately 25 m the probability of mortality within analysed drought period 2015-2019 exceeds 50% and grows significantly up to the highest observed SI of 40 m (Fig. 7d).
The simultaneous effect of stand age and SI on the probability of mortality is illustrated by a matrix showing mortality probability in the period 2015-2019 for 10-year age classes of stands and 1-metre SI ranges (Fig. 8). For stands up to 50 years of age, the probability of mortality was close to 0 independent of SI. In stands of 60 to 80 years of age, a relatively small risk of mortality (2-27%) occurred only in stands growing in the most productive sites (SI> 25 m). All stands over 100 years of age were at risk of mortality but on less productive sites (SI up to 25 m), the probability of mortality was up to about 30%, while for the most productive sites (SI over 30 m), the probability of mortality exceeded 50% for most of the stands.
The model with incorporated spatial coordinates explained 47.1% deviance and was used for developing the map of probability of mortality for all Scots pine dominated stands within the study area (Fig. 9). The parameters of the model incorporating spatial coordinates are presented in the Supplementary Table 2.

Discussion
Our results show that there is a one-year lag between drought occurrence and mortality. Furthermore, we have shown that, despite the multifaceted stochastic nature of tree mortality, during the drought the probability of mortality is primarily determined by stand age and the site productivity. Paradoxically, forests growing on the poorest, including the driest, sites with SI up to 20 m, are comparably less vulnerable to mortality. Greater mortality probability is observed in more productive stands and increases with SI. The probability of mortality is also increased by stand density and modified by aspect.
A reasonable explanation for the delayed drought effect is the sequential phenomenon of droughts, which is recognised as a trigger for mortality in the boreal forest (Sánchez-Pinillos et al., 2021). In a drought year, the tree is weakened and the following drought year may lead to more frequent mortality. This may explain why mortality was not as severe in 2015, one of the most intense drought years, and increased in subsequent drought years. Our results indicate that increasing stand age in interaction with higher site productivity resulting simultaneously in increased stand height leads to enhanced susceptibility to mortality observed in drought period 2015-2019. Amongst the stands with mortality, there were those in which no pests were observed or mainly two pests were documented: the sharp-dentated bark beetle (Ips acuminatus) and the fungus Sphaeropsis sapinea. The unambiguous identification of drought or pests as the direct cause of mortality is limited (Mezei et al., 2017). However, the indications presented in our study suggest that drought was the main cause of mortality. Firstly, as shown in the first stage of the analyses, the scale of mortality was closely related to CWB. Secondly, the pest species contributing to mortality belong to the group of secondary pests (Plewa and Mokrzycki, 2017;Skrzecz andPerlinska, 2018, Sikorski et al., 2023). Their gradation occurs after the stands have been previously weakened by drought (Jactel et al., 2012;Rouault et al., 2006). The fact that the scale of mortality in a given year was dependent on the CWB in the previous year further indicates the relationship with  stand weakening and the secondary nature of pest gradation. Moreover, prior to the sequence of years with dry and hot summers in 2015-2019, both of these pests occurred very sporadically and on a small scale (Jankowiak et al., 2021;Jankowiak and Bilański, 2013;Kwaśna et al., 2017;Plewa and Mokrzycki, 2017;Sukovata et al., 2021). Although CWB was also slightly negative for the period 2009-2014, no significant correlation of mortality with stand age and site productivity was observed for this period. This may be due to the different severity of the drought expressed by CWB (Fig. 4). At lower drought severity, the stands with largest and oldest trees, which tend to have deeper root systems, may be relatively drought tolerant compared to younger and smaller trees (Brunner et al., 2015;Padilla and Pugnaire, 2007). These results are consistent with recent research on the physiological tree responses to drought (Adams et al., 2017). Drought-induced tree mortality is associated with two physiological mechanisms: (1) loss of xylem function due to embolism, which partially or completely reduces water transport through the conductive system and leads to tissue desiccation (Choat et al., 2018;Olson et al., 2018); (2) carbon starvation involving a loss of balance between non-structural carbohydrate requirements and supply, which limits the ability to meet tree requirements related to metabolic, osmotic, and defense functions (Adams et al., 2017;Sevanto, 2018). We show that the oldest and highest Scots pine stands growing on the most productive sites exhibit increased mortality susceptibility. We argue that this is biologically sound for two reasons. Firstly, water-conducting conduits of the highest trees from the most productive sites are more vulnerable to xylem embolism (Adams et al., 2017;Olson et al., 2018). Moreover, taller plants are also more susceptible to carbon starvation (Adams et al., 2017;Choat et al., 2018), contributing to the weakening of trees, which may result in increased mortality with some delay. The ensuing depletion of non-structural carbohydrate reserves leading to mortality is specific to taxonomic groups. In gymnosperms, carbon starvation is closely connected with xylem hydraulic vulnerability and of far greater importance than in angiosperms (Adams et al., 2017). Trees that grow faster optimise their anatomical structure for more efficient water transport while increasing allocation to branches at the expense of fine roots (Schuldt et al., 2016). This optimisation is diametral to physiological optimisation to resist drought stress (Jacobs et al., 2021). Also, assimilate transport is strongly connected with water transport (Adams et al., 2017;Hesse et al., 2021;Sevanto, 2018). Consequently, especially in large mature trees due to the need to transport assimilates over long distances, imbalances of non-structural carbohydrates in fine roots may occur during drought (Landhäusser and Lieffers, 2012;Ryan et al., 2014). Moreover, large, fast-growing trees on most productive sites show increased water demand resulting from relatively higher metabolic requirements (Aranda et al., 2012;Meinzer et al., 2009). Consequently, mature trees growing faster may be more susceptible to increased mortality (Jacobs et al., 2021).
Secondly, the ability of trees to adapt to changing environmental conditions, including the plasticity of the root systems, may also play a significant role (Brunner et al., 2015). The ability of the tree root system to adapt to site conditions decreases with age. The taproots of Scots   Fig. 6. Partial effects of individual explanatory variables (stand age, site index, stand density and aspect) on the drought-induced stand-level mortality probability. The green area indicates 95% confidence interval for the spline function of GAM model (red line). Probability refers to the risk of mortality in a given stand covering at least 0.1 ha.
pines reach their maximum depths after approximately 30-50 years (Correa et al., 2019). Therefore, only younger stands can fully adapt their root systems to reductions in soil moisture. Moreover, elongation of vegetation period resulting from climate change, nitrogen deposition and increased CO 2 concentration accelerate aboveground biomass growth and proportionally reduce the size and depth of the root system, which may further increase susceptibility to drought (Brunner et al., 2015;Schuldt et al., 2016). Apart from increased probability of mortality during drought, decreased growth resistance and reduced resilience are also observed in older stands. Based on a meta-analysis from 23  Site index is an indicator of site productivity that was determined based on stand height calculated for a base age of 100 years. Boxes for which there were no stands with a given value of the age -site index combination are marked in white. The probability shown in the graph represents the average probability that a stand in a given age range and site index experienced mortality in 2015-2019 on an area of at least 0.1 ha. experiments, Sohn et al. (2016) showed decreased growth resistance with increasing stand age during drought. And specifically, for Scots pine, Martínez-Vilalta et al. (2012) showed reduced resilience to drought of older trees. However, the mortality patterns are related to a complex set of factors, and, at the same time, they are species-or taxa-specific (Stephenson and Das, 2020). Also because of the managed nature of central European forests, data from real old-growth forests is missing in this and similar analyses.
We also found a positive correlation between mortality probability and stand density. Lower stand density alleviates the adverse effects of drought by reducing competition between individual trees (Andrews et al., 2020) and improves the water balance, reducing overall water stress. Drought vulnerability in less dense managed pine stands is lower than in unmanaged stands, experiencing higher competition for soil water (Lucas-Borja et al., 2021). Research by Sohn et al. (2016) andD'Amato et al. (2013), showed that reducing stand density by thinning interventions improves the adaptation of pines to drought and recovery of growth after drought in young forests. The mortality probability may also be modified by the structure and species composition of forests. We analysed stands with a predominance of Scots pine, primarily even-aged with a single-storey structure. Different stand density effects can be expected for complex, mixed and uneven-aged stands (Pardos et al., 2021).
The modifying effect of the aspect on the probability of mortality is also most likely caused by the different soil moisture resulting from the differences in temperature distribution and radiation (Iverson et al., 1997). Also Bergh et al. (2005) found that the strongest periodic water shortages occur on south-western facing slopes, which explains the increased mortality of stands located in the azimuth range from 200 to 250 (Fig. 6). Moreover, we found a significant increase in explained deviance -up to 47.1% -when coordinates were included in the model. An explanation for the large coordinate influence may be that drought-driven forest mortality is patchy and stands are more susceptible in some areas than in others (Tague, 2022). This is most likely due to the influence of factors not included in the model, which may be related to both stand and site characteristics, groundwater level, management history and provenance Callahan et al., 2022;Soong et al., 2020).
The spread of remote sensing technologies allows continuous and precise mapping of forest ecosystem attributes over large areas (Jurjević et al., 2020;Tymińska-Czabańska et al., 2022). Thanks to the full coverage of the study area with ALS data, it was possible to characterise the individual attributes of tens of thousands of forest stands. Obtaining such an extensive data set is impossible using data from sample plots. ALS data are particularly useful for determining top height and the number of trees in stands (White et al., 2013;Wang et al., 2019;Jurjević et al., 2020). Certainly top height of stands determined from sufficiently dense ALS data are much more accurate than heights measured by traditional field inventory (Wang et al., 2019). When determining stand density on the basis of ALS, the problem is to calculate the number of trees in complex, especially broadleaved stands and the number of trees in the lower layers of the stand (Socha et al., 2020a;Iqbal et al., 2021). However, in the case of the Scots pine monocultures we studied, the number of trees determined by ALS can be considered quite precise. To avoid errors associated with underestimating the number of trees of the lower storey, we characterised stand density using an approach that considered only trees from the upper layer of the stand according to the approach proposed by Socha et al. (2020a).
Our results extend and add new insight to earlier studies linking increased mortality and larger tree sizes (Brienen et al., 2020;Olson et al., 2018;Stovall et al., 2019) by documenting a relationship between increased mortality and stand age, site productivity and stand density. Our research indicates that younger stands are less exposed to drought-induced mortality than older stands. Moreover, to the best of our knowledge, previous studies have not documented a direct link between site productivity and increased mortality. Yet, this link reveals a new dimension of the widely observed forest site productivity increases worldwide (Albert and Schmidt, 2010;Appiah Mensah et al., 2021;Metslaid et al., 2011;Pretzsch et al., 2014b;Socha et al., 2021b) and could help to link these observations to seemingly contradictory findings of increasing forest mortality . Thus, in light of our results, increasing site productivity is not necessarily a positive phenomenon enhancing carbon uptake and forest yields but might as well have serious ecological consequences for forest ecosystem functioning and service provisioning. As a result of significantly increased site productivity in recent decades, stands have reached larger heights and densities (Pretzsch et al., 2014a(Pretzsch et al., , 2014b. Periods of favourable growing conditions that accelerated tree growth may have led to a structural overshoot of aboveground tree biomass, which, under drought conditions, results in a subsequent temporal mismatch between water demand and water availability (Jump et al., 2017). This results in increased water stress, which contributes to accelerated tree mortality (Jump et al., 2017). Elevated stand-level mortality connected with SI and density means that for a given drought intensity, forests are now at higher mortality risk than they were a few decades ago. Distinguishing between the effects of drought and pests on stand mortality is very difficult due to the strong link between pest occurrence and drought. Therefore, in forecasting mortality risk for forest management made for different climate scenarios, it seems reasonable to treat these causes of mortality together. However, it must be taken into account that long-term forecasting of forest mortality is nevertheless subject to great uncertainty due to the very high variability of this phenomenon. An additional source of uncertainty may be that, in addition to the climatic conditions, particularly climatic water balance, the risk of mortality is linked to changes in site productivity that are difficult to predict in the long term.
To avoid disturbing the results by taking into account the selfthinning of individual trees, as stands with mortality we defined those in which the area of forest cover loss due to sanitation felling was larger than 0.1 ha. The adoption of the solution with a minimum area instead of a minimum volume was justified by the fact that if the same proportion of trees dies in terms of volume or number, the probability of mortality in all stands, regardless of their volume and density, will be the same. Potential uncertainties arising from the use of sanitation felling data are that, in addition to dead trees, trees showing strong symptoms of mortality but still alive may also be harvested, decreasing the Fig. 9. Probability of mortality as a function of stand age, site index, stand density, aspect, CWB and coordinates predicted for all Scots pine dominated stand within the study area.
accuracy of estimation of the volume of dead trees. However, this limitation has no significant effect on the validity of the determination of the presence of mortality in a given stand. An additional limitation may be that individual trees that died in a given year may have been missed by the forest monitoring and removed the following year. However, on the other hand, it is highly unlikely that such cases occurred for the areas of at least 0.1 ha that we used as a threshold in the analyses. Another limitation is that our approach allows the identification of stands with higher mortality risk but does not allow a precise calculation of the volume, share and number of dead trees. Further research is therefore needed to quantitatively describe mortality risk at the level of individual stands, which, despite the development of drought research, including physiological studies, is still difficult (Trugman et al., 2021). Our results refer to the stand level. To draw conclusions about individual trees, further research utilising local data at the individual tree level are needed. Global and regional vegetation models predict an increase in vegetation carbon in the future (Friend et al., 2014) yet they struggle to properly account for mortality processes and the crucial role of carbon residence time (Bugmann et al., 2019;Friend et al., 2014). Uncertainties in carbon residence times might result from increasing frequency of disturbances (Senf and Seidl, 2020) but we show that they might also result from increased mortality susceptibility due to the interaction of improved growing conditions and drought stress for older and most productive forests. Therefore, the risks associated with forest growth stimulation by climate change resulting in a lagged increase in canopy tree mortality and reduced carbon storage by forest ecosystems should be considered when modelling vegetation carbon pools (Brienen et al., 2020;Büntgen et al., 2019;Körner, 2017;Xu et al., 2021).

Conclusions
We found that the occurrence of mortality in Scots pine stands is mainly driven by the lagged effect of the water deficit in the period May-August of the previous year expressed by the CWB. Moreover we documented that the oldest stands growing on the most productive sites have the highest susceptibility to increased mortality during drought. Climate change increases the frequency and severity of drought and, however, simultaneously along with nitrogen deposition, increases site productivity and accelerates forest growth. Furthermore, the predicted increase in temperature associated with climate change will therefore lead to a decrease in CWB, even with similar levels of precipitation. Therefore, increasing site productivity and decreasing CWB may accelerate the intensity and prevalence of drought-induced forest mortality. Exacerbated susceptibility to mortality should be considered in forest carbon sink projections, forest management, and policies designed to increase resilience and protect forest ecosystems. Our results are important to better understand forest ecosystem functioning and ecosystem service provisioning under climate change and to plan the adaptation and the sustainable forest management. Several recent studies indicate that the increase in global forest mortality is one of the most pressing issues today. Therefore a deeper understanding of the drivers and factors underlying forest mortality is urgently required to identify management options and understand implications for forest ecosystem service provisioning such as timber supply, carbon sequestration and providing habitat for biodiversity.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The authors do not have permission to share data.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.agrformet.2023.109680.