Stand characteristics modulate secondary growth responses to drought and gross primary production in Pinus halepensis afforestation

Impacts of climate warming on forests vigour are forecasted to increase in magnitude. Yet it remains unclear how stand characteristics and competition modulate the relationship between tree growth and gross primary production with drought. Here, we studied how the spatial variation in stand density, basal area and height modulates tree growth (Basal Area Increment, BAI and stand growth), summer NDVI, as well as their responses to drought (Standardized Precipitation-Evapotranspiration Index, SPEI) in 56 Aleppo pine (Pinus halepensis Mill.) planted forests located in Northeast Spain. Long-term BAI responses to SPEI were strongly determined by stand density, suggesting that competition modulates Aleppo pine growth responses to drought. Along this, summer NDVI also displayed strong associations with SPEI. NDVI was mostly related with stand growth, suggesting canopy densification drives NDVI pattern and trends. Short-term BAI and NDVI responses to severe droughts were mainly independent of stand characteristics. In the studied region, drought is a universal factor limiting Aleppo pine secondary growth and canopy greening. However, the results suggest that stand density modulates Aleppo pine growth responses to drought on the long-term, reducing the growth in densest stands. Denser stands with larger trees are the ones that present higher NDVI values, suggesting that canopy activity depends more on stand canopy coverage than on secondary growth rate and its response to drought. In these Mediterranean pines, canopy activity and secondary growth are temporally coupled but spatially decoupled.


Introduction
As climate becomes warmer, the likelihood of hotter drought occurrence and their impacts on forest vigour increases (Allen et al. 2015). This trend may be particularly dramatic in the Mediterranean region, where warming rate is far higher than average, threatening forest persistence in this drought-prone area (Cramer et al. 2018). Thus, monitoring forest growth by combining different sources of information becomes paramount to anticipate Mediterranean forest response to climate change (Gao et al. 2020;Ruiz-Benito et al. 2020).
Individual tree basal area increment is a reliable proxy of tree vigour, enabling to understand the response of forests to climate change (Xu et al. 2017;Marchand et al. 2019). Treering series are valuable tools to monitor both growth trends and tree responses to drought (e.g. Sánchez-Salguero et al. 2010;Camarero et al. 2015;Babst et al. 2018;DeSoto et al. 2020). Both, climate-growth relationships, and the capacity of trees to resist severe droughts and to recover after them (sensu Lloret et al. 2011), inform on how trees respond to drought at long-and short-term scales, respectively. Variations in stand characteristics can modulate the trees response to drought. For example, tree growth is age-dependent (Martínez del Castillo et al. 2018) and so growth trends and their Communicated by Lauri Mehtätalo. 1 3 response to drought can vary along time (Lloret et al. 2011) with trees being more vulnerable to drought at small (Colangelo et al. 2017) and large sizes (Olson et al. 2018). Competition between neighbouring individuals has been found to have profound impacts on forest growth and its response to climate (Gea-Izquierdo et al. 2009;Linares et al. 2010;Rozas 2015). Along this, dense stands have been found to be more vulnerable to drought (Bottero et al. 2017;Tsamir et al. 2019;Andrews et al. 2020) and thinning is a wellknown management strategy to increase stand growth resilience to drought (Olivar et al. 2014;Manrique-Alba et al. 2020). In dense stands, negative interactions modify canopy architecture and light interception and can ultimately lead to tree mortality (Olano et al. 2009). Stand characteristics and tree-to-tree interactions can modulate secondary growth responses to drought, which can be quantified by measuring and analysing tree-ring growth patterns. However, upscaling ring-width measures to understand trends in forest productivity is challenging Kannemberg et al. 2019;Vicente-Serrano et al. 2020b).
Remote-sensing data products offer information on how canopy greening, and thus gross primary production respond to drought and other factors along time (e.g. Senf et al. 2018). The Normalized Difference Vegetation Index (NDVI) is a measure of the photosynthetically active radiation obtained by vegetation and thus it relates to the Gross Primary Production (GPP) (Myneni et al. 1995). Tree ring information has been correlated to NDVI (e.g. Vicente-Serrano et al. 2020a), enabling our ability to upscale secondary growth patterns at regional scale and providing both estimates of carbon gain and information about forest resilience to drought Gazol et al. 2018). However, secondary growth can also be decoupled from GPP if other C sinks, such as survival, reproduction, or defence, are prioritized after drought occurrence (Doughty et al. 2015;Kannenberg et al. 2019).
The temporal scale at which NDVI and radial growth maximize their relationship and their spatial covariation have received substantial attention over the last decades (e.g. Vicente-Serrano et al. 2016, 2020bXu et al. 2019). In this sense, Xu et al. (2019) found a strong positive relationship between summer NDVI and ring-width index indicating that both, canopy greening and tree growth disruptions and recovery patterns were related to water availability. NDVI and radial growth relationship is contingent to scale phenomena, i.e. the resolution of remote sensing data products (i.e. pixel size) as well as the representativeness of tree-ring sampling strategies (Nehrbass-Ahles et al. 2014;Babst et al. 2018). However, little attention has been shed on how local factors modulate the relationship between secondary growth and NDVI with climate (Vicente-Serrano et al. 2020b).
Pine plantations in the Iberian Peninsula offer a valuable opportunity to study the relationship between radial growth and NDVI responses to drought in the Mediterranean region, as they are monospecific and even-aged, and stand structure can be estimated accurately through remotely sensed data (Aragones et al. 2019). This is the case of the large number of Aleppo pine (Pinus halepensis Mill.) plantations established in the eastern part of the Iberian Peninsula during the second half of the past century (Serrada 1992;Montero 1999). This species was favoured because it displays several strategies to tolerate drought, such as a plastic response to water availability (De Luis et al. 2009;Camarero et al. 2010). Intra-annual patterns of stem radial increments reveal a bimodal growth pattern of the species with peaks in spring and autumn and a rest during the dry summer (Pacheco et al. 2017). The bimodal growth pattern points to the importance of favourable climatic conditions during winter and early spring to the carbon gain of the species (Maseyk et al. 2008;Camarero et al. 2010). Dendrochronological studies also evidenced that wet conditions during winter and spring favour the secondary growth of Aleppo pine (Pasho et al. 2011), whereas severe drought events may result in widespread growth declines (e.g. Gazol et al. 2018;Helluy et al. 2020). Despite this species ability to recover pre-drought growth rates after severe droughts , there are also evidences of drought-induced mortality episodes in planted stands of dry regions or in populations located under xeric conditions (Sánchez-Salguero et al. 2010;Dorman et al. 2015). Analyses of NDVI time series also indicate that Aleppo pine canopy activity decreases as spring advances and stops during the summer months (Aragones et al. 2019). However, it remains unclear how changes in secondary growth responses to drought translate into NDVI variations, and to what extent stand characteristics modulate these relationships in Aleppo pine planted forests.
In this study, we evaluate how the stand structure affects secondary growth and NDVI response to drought in a network of semi-arid Aleppo pine afforestation. To this end, we use a robust dataset comprising annually resolved ring-width series of 56 Pinus halepensis stands in Northeast Spain and field measurements of stand characteristics: stand density, height, diameter and site quality. Ring-width measures were used to quantify basal area increments (BAI) in each stand for the period 1986-2016. Further, we used medium-resolution NDVI data (30×30 m) to determine canopy activity and gross primary production variations in the period 1986-2016. The Standardized Precipitation-Evapotranspiration Index (SPEI) was used to characterize drought severity . We evaluated how BAI and NDVI were related with SPEI, and how severe drought events affected the resistance and recovery capacity of both BAI and NDVI. Besides, we considered how stand characteristics modulate these relationships. We hypothesize that (i) there are significant relationships between BAI and NDVI with SPEI, reflecting that both, secondary growth and canopy greening respond to drought; (ii) these relationships are mediated by stand characteristics particularly because neighbouring competition can limit secondary growth (BAI) and modulate canopy architecture (NDVI); and (iii) severe drought events cause abrupt BAI reductions and drops in NDVI, but stand characteristics modulate how the stands respond to and recover after them. In the complex process of tree growth, several steps are involved from canopy development to secondary growth dynamics, and thus NDVI and tree-ring width measures can provide complementary information on forest response to drought.

Study sites and species
The study site comprises 56 Aleppo pine plantations located in the Ebro basin (Aragón autonomous region, Northeast Spain) ( Fig. 1), which is characterized by semi-arid Mediterranean continental conditions. Mean annual temperature oscillates between 11 and 16 ºC and total annual precipitation oscillates between 300 mm in the valley bottom and over 800 mm in the surrounding mountain areas (Vicente-Serrano and Beguería 2003). The driest period occurs during summer, and long dry periods are frequent due to the high potential evapotranspiration and a high interannual variability in precipitation (Vicente-Serrano and Beguería 2003;Vicente-Serrano et al. 2014).
Stand density (trees·ha −1 ) varied considerably from near 200 to 3000 individuals per hectare as well as stand basal area (Table S1). Similarly, site index, which indicates the quality of the stand in terms of productivity, differed between sites. Despite the presence of some young and old stands, most sites presented ages ranging between 50 and 60 years (Table S1).

Climate and NDVI data
We used 1.1 km 2 -gridded climatic data for precipitation and potential evapotranspiration assembled at a weekly timescale based on values of daily mean temperature, sunshine duration, wind speed, relative humidity, and precipitation coming from the Spanish Meteorological Agency (AEMET). Time series for each stand in the period 1961-2018 were gattered based on the location of the stands within the gridded climate data. The dataset was homogenized and quality-checked prior to the estimation of gridded values (see Vicente-Serrano et al. 2017). The authors estimated potential evapotranspiration according to the FAO-56 Penman-Monteith (PM) equation which is a widely used method (see Vicente-Serrano et al. 2017). Monthly precipitation (P) and potential evapotranspiration (PET) values were aggregated by summing weekly values in the gridded dataset. Finally, we calculated the monthly water balance (P-PET) as the difference between precipitation and potential evapotranspiration. Further, the water balance was used to calculate the Standardized Precipitation-Evapotranspiration index (SPEI), which is a drought index that allows for spatial and temporal comparability of drought severity (Vicente-Serrano et al. 2010). SPEI is a multiscalar drought index that can take negative and positive values, indicating dry and wet periods, respectively (Vicente-Serrano et al. 2017). The temporal scale is determined by the number of months considered to calculate the water balance. As an example, the 3-month long SPEI for March takes into consideration the water balance in January, February and March. We calculated 3, 6, 9 and 12-month long SPEI for March, June, September and December in each stand for the entire study period . To calculate the SPEI we used the SPEI R library (Beguería et al. 2014).
The Normalized Difference Vegetation Index (NDVI) data were obtained from the Landsat ETM + , TM, and OLI image collections held in Google Earth Engine (Gorelick et al. 2017). We generated 56 NDVI dense time series for the period 1986-2016, one per study stand, by averaging pixel (spatial resolution of 30×30 m) values around the field plot centre. Each time series has a temporal resolution 16-8 days depending on epoch and satellite availability. Observations with atmospheric contamination (e.g. clouds, shadows) were vetted via the Fmask algorithm (Zhu and Woodcok 2015). For each year, we calculated the average NDVI value in each stand during the spring (March to May) and summer months (June to August). The two values were strongly correlated (r = 0.80; p < 0.01) and so we focused on summer months which correspond to the period with the lowest values across years (Fig. S1).

Field sampling and dendrochronological data processing
Fifty-six circular plots of an area of 225, 400 or 625 m 2 , depending on tree density, were established to characterize the stands. The number of adult trees (diameter at breast height, dbh ≥ 7.5 cm) was counted and their dbh measured to estimate stand density (n ha −1 ) and basal area (m 2 ha −1 ), respectively. An increment core from the bark to the pith of 5 to 6 individuals was extracted at 1.3 m height with a Pressler increment borer (5 mm; Haglöf, Sweden) for further dendrochronological processing, leading to a total of 333 trees. The height (m) of the selected trees and, if needed, of an ancillary sample of trees to estimate Assmann dominant height, was quantified using a Haglöf Sweden®, Vertex instrument. Moreover, an increment core from the bark to the pith was extracted at the base of two dominant trees to estimate stand age. Table S1 shows the main summary statistics of the stand variables for sampled plots, and Table S2 factors affecting their site-to-site variation.
Stand structure was quantified using stand density, stand basal area and site index. Stand density was measured as the number of sampled trees per hectare by counting the number of trees in each plot multiplied by the area sampled. Stand basal area (m 2 ha −1 ) was measured as the basal area of the sampled trees divided by the number of trees present in the stand and extended to hectare. Stand density and stand basal area (m 2 ha −1 ) were considered to account for differences in the degree of inter-tree competition between stands (Wright et al. 2018). Both variables vary between sites but not in time as they were measured only once. This approach misses the contribution of the trees that have died prior to our sampling, but it conveys a reasonable estimate of between stand differences. Since tree age can modify growth trajectories substantially (Martínez del Castillo et al. 2018), we also considered potential effects of estimated tree age.
Additionally, we used site index as a measure of site quality. To do that we followed Rojo-Alboreca et al. (2017) which constructed polymorphic, base-age-independent site curves based on dominant tree height. Estimated stand age was quantified as the maximum age of the trees sampled in each stand. The site index informs on the quality of the stand according to dominant tree height (Rojo-Alboreca et al. 2017).
The extracted cores were air-dried, glued and polished using a series of sand-paper grids until tree-ring boundaries were clearly visible. After that, they were visually crossdated using standard dendrochronological methods and measured under magnification to the nearest 0.01 mm using a measuring stage (Velmex Inc., Bloomfield NY, USA). The COFECHA software was used to check cross-dating accuracy by calculating moving correlations between each treering series and the mean study site series (Holmes 1983).
In total 333 Aleppo pine trees were accurately cross-dated. The mean inter-series correlation was 0.212 and the average mean sensitivity (Cook and Holmes 1996) of the dated trees was 0.181.
Tree-ring width (mm) measures were transformed into basal area increments (cm 2 year −1 ) according to Biondi and Qeadan (2008). Basal Area Increment (hereafter BAI) measures were averaged year-by-year across individuals within stands and multiplied by stand density to estimate standlevel growth (hereafter stand growth). Thus, for each stand we measured BAI (cm 2 year −1 ) at the individual tree level and estimated stand growth (m 2 ha −1 year −1 ).
Short-term responses to drought were quantified as the BAI and summer NDVI reduction observed in years with abnormally low values of SPEI in the period 1986-2016. Particularly, we selected the years 1994, 2005 and 2012 because they presented the lowest values of 12-month SPEI (< − 1.5 in all cases) in the period (see Fig. 2). Drought impacts were quantified as drought resistance (Rt; sensu Lloret et al. 2011). The 'resistance' is conceptually identical to the concept of abrupt growth changes described in Schweingruber et al. (1990) and was calculated as the ratio between the growth in the drought years and the growth in the preceding n-years. In this study, we used a period of three years before the drought to calculate the resistance. To quantify the capacity of BAI and summer NDVI to recover after drought we focused on two indices: the recovery period and the relative average growth reduction (Thurm et al. 2016;Schwarz et al. 2020). The recovery period quantifies the time required (in years with decimals) to reach pre-drought levels of growth (Thurm et al. 2016), whereas the relative average growth reduction quantifies the cumulative growth reduction during the recovery period divided by the length of the recovery period (Schwarz et al. 2020). Average growth reduction was expressed in percentage (i.e. relative average growth reduction). We considered a period of 4 years as the maximum length of the recovery period because that was the time between the first year after the last drought considered (2013) and the end of the BAI and NDVI studied periods (2016). In the case of BAI, calculations were done at the tree level while in the case of NDVI calculations were done at the stand level.

Statistical analyses
Our main aim was to evaluate whether the stand structure modulated the relationships between BAI, stand growth and NDVI with drought (SPEI). For this purpose, we used a series of Linear-Mixed Effect Models (LME) (Pinheiro and Bates 2000) using BAI, stand growth and NDVI for the period 1986-2016 as response variables. In general, models were of the form: where Y is the response variable (either BAI, stand growth or NDVI), f(x) is the set of fixed effects, u s represents the plot (or tree nested within plot) normally distributed random effects, v t is a normally distributed random effect for calendar year t (year effect), and e st is the normally distributed residual for tree (stand) s at year t. We used this model adapted from the one proposed by Mehtätalo et al. (2014) to study the temporal variation in tree growth between stands while accounting for the effect of unspecified tree-, plot-and year-level factors. In the case of BAI, we used nested tree effects within plots as a random factor as tree ring-width series represents repeated observations over the same individuals in each stand. In the case of stand growth and NDVI, we used plot as random factor. In all cases, year was used as a crossed random factor to account for unmeasured temporal variations.
First, we created a model including stand growth as response variable and NDVI as fixed effect to study the relationship of these variables in the period 1986-2016. In this case, stand and year were used as random factors (random intercept) and stand growth was log-transformed (log(x + 1)).
Second, we studied the relationship between BAI and NDVI with climate and stand characteristics. The fixed part of the model consists of the age trend (only in the case of BAI), the drought index and stand characteristics: stand density, stand basal area and site index. Preliminary analyses showed that BAI had strong positive relationships with SPEI calculated at different temporal scales but the strongest relationship was found when the SPEI was calculated at 12-month scale for the month of September. In the case of BAI, the age trend was represented by tree age as fixed factor in the models. Both, BAI and tree age were log-transformed prior to the analyses. We included potential interactions between SPEI with stand density, stand basal area and site index. We prevent models containing both, stand density and stand basal area as these variables are related. Each model was evaluated graphically, and its fit was quantified by computing the marginal (due to fixed effects) and conditional (due to fixed and random effects) pseudo-R 2 (Nakagawa et al. 2017). Model selection was performed by comparison of the Bayesian Information Criterion (BIC) of the different models (Burnham and Anderson 2002) using ML and best model parameters were recalculated with REML (Zuur et al. 2009).
We also studied the temporal variations of the short-term response of BAI and NDVI to three selected severe droughts (1994, 2005 and 2012), and whether stand characteristics modulate stands response to severe droughts as well as postdrought recovery. Short-term responses of BAI and NDVI to drought were quantified as the resistance, recovery period Y = f X st + u s + v t + e st and the relative average growth reduction. For each variable, we created LME models that contained a factor differentiating the three droughts (i.e. 3-levels, 1994, 2005 and 2012, with SPEI value < − 1.5) as fixed factor. Stand identity, in the case of NDVI, and tree nested within site in the case of BAI, were used as a random factor, and response variables were log-transformed (log(x + 1)) if required. As before, we included potential interactions between the SPEI with site index and stand characteristics to decipher if site-to-site variations in response to SPEI were different between droughts. Differences in short-term responses of BAI and NDVI to drought between droughts (1994, 2005 and 2012) were estimated as Least-square means calculated for the final model (Lenth 2022). All statistical analyses were performed in R environment (v. 3.5.2 The R Foundation for Statistical Computing 2018). The "dplR" package (Bunn et al. 2018) was used to calculate tree growth (BAI), the "pointRes" was used to calculate resilience components (van der Maaten-Theunissen et al. 2015) and the "lme4" (Bates et al. 2015) and "MuMin" (Barton 2012) packages were used to calculate the linearmixed effect models and perform the multi-model selection, respectively. The "emmeans" package was used to calculate least-square means (Lenth 2022) and the "effects" package was used to visualize models (Fox 2003;Fox and Weisberg 2019).

Results
Aleppo pine secondary growth (stand growth) presented a temporal pattern tightly coupled with the temporal variations in SPEI (Fig. 2) as well as the temporal trend in NDVI (Fig. 2b). NDVI also showed a positive trend (Fig. S2). The temporal variation in stand growth and summer NDVI were positively related (β = 0.217 ± 0.026; t = 8.289; R m 2 = 0.089; R c 2 = 0.607; Fig. 3). The model relating individual tree BAI with drought (12-months long September SPEI) and stand characteristics (stand density, stand basal area and site index) revealed the simultaneous effect of both factors (Table 1; Table S3). Site index had a positive impact on BAI while stand density had negative impacts (Table 1). Notably, the effect of SPEI was dependent on stand density (Table 1) suggesting the existence of more plastic growth responses to SPEI in sites with lower density (Fig. 4). Interestingly, we found a clear positive age effect indicating clear growth trends.
The results of the models testing for the relationship between NDVI with drought (12-months long September SPEI) while accounting for stand characteristics pointed to the importance of stand basal area as modulator of canopy greenness (Table 1; Table S3; Fig. 5). Stand basal area has a significant positive impact on NDVI (Table 1; Fig. 5).
Similarly, NDVI increased in wet years as suggested by the positive associated with SPEI (Table 1).
Short-term BAI resistance to drought and their capacity to recover after a dry period was variable (Fig. 6) but did not differ depending on stand structure (Table 2 and Tables  S4-S5). The resistance facing drought in terms of BAI was lower (Table 2) in the 2005 than in the 1994 and 2012 droughts. Consistently, the recovery period and the relative average growth reductions after the 2005 drought were higher than after the 1994 and the 2012 droughts ( Fig. 6; Table S5). There were 43, 170 and 15 trees in which the BAI did not recover to pre-drought growth levels in the four years following the 1994, 2005 and 2012 droughts, respectively. Contrary to BAI, the resistance facing drought in terms of NDVI was lower in the 1994 and 2012 droughts than in the 2005 drought (Fig. 6). Consistently, the relative average NDVI reductions after the 2005 drought were lower than after the 1994 and the 2012 droughts (Fig. 6). The recovery period after the 2012 drought was significantly lower than the recovery period after the 1994 and 2005 droughts. Short-term NDVI responses to drought were independent of stand characteristics (excluding relative growth reductions, Table 2 and Table S4, S5) but they varied depending on the variation in SPEI between sites (Table S5).

Discussion
The results provided in this study indicate that Aleppo pine secondary growth and canopy greening are strongly sensitive to water availability, corroborating our first hypothesis.  Besides, we found an interaction between stand density and SPEI on driving BAI, suggesting that long-term growth responses to SPEI are associated to stand density, being more plastic in less dense stands. These results align with the already known sensitivity of Aleppo pine growth to drought and suggest that secondary growth responses to climate are mediated by neighbouring competition (Linares et al. 2010;Rozas 2015). Our results also show that NDVI presents higher values in stands with higher basal area, and that it has increased over the period 1986-2016 in the studied stands. Neither tree resistance to drought, nor the recovery capacity, were related to stand variables, which lead us to reject our third hypothesis. All together, these results suggest that regulations of stand density can increase the growth capacity of planted Aleppo pine stand in normal years but has no clear effects on modulating the capacity to withstand rare and severe droughts. That is, when drought intensity surpasses certain threshold its impacts on growth cannot be ameliorated by solely modifying stand density. The question that emerges here is whether potential mortality consequences caused by severe droughts could have been lower in less dense stands. Aleppo pine growth is sensitive to drought (De Del Castillo et al. 2015). In our study, this is supported on: (i) the strong linkages between long-term growth trends and SPEI (Figs. 2, 4); and (ii) the abrupt short-term growth reductions observed during the occurrence of severe droughts in the years 1994, 2005 and 2012, particularly in 2005 (Fig. 6). Pasho et al. (2011) found that Aleppo pine growth responds to drought more tightly at scales of 6 to 15 months. (i.e. including both the current year and the year prior to ring formation) indicating that tree growth is better explained by SPEI at yearly scales. In accordance, we found a strong impact of the SPEI calculated at scales of 12-months in September (Fig. 4), pointing to the importance of soil water reserves Del Río et al. 2014). We also found that summer NDVI responded to the SPEI suggesting that soil water reserves play an important role for canopy greening. Thus, synchrony between secondary growth and canopy activity partially depend on growth response to drought (Xu et al. 2019) but can also be affected Fig. 4 The secondary growth of the sites a is separated according to the stand density (low-below the 30 percentile; medium-sites between the 30 and the 70 percentile; high-above the 70 percentile; Table S1), as well as the effect of the interaction between the SPEI at the 12-month scale for September and stand density on BAI according to the models proposed in Table 1 b The brown, light blue and dark blue colours correspond to low, medium, and high density, respectively

Fig. 5
Relationship between summer NDVI and the SPEI at the 12-month scale for September. The solid line represents the effect according to the linear mixed-effect model (Table 1) whereas the dashed lines indicate the 95% confidence interval for the effect. Sites are separated according to their stand basal area: low-below the 0.3 percentile; medium-sites between the 0.3 and the 0.7 percentile; high-above the 0.7 percentile by variations between the local conditions of the stands not considered in this study.
Secondary tree growth was greater in stands with lower density suggesting a potential influence of neighbouringcompetition on growth (Linares and Camarero 2010;Rozas et al. 2015;Andrews et al. 2020), but also in sites of higher site quality, pointing to the positive effect of edaphoclimatic conditions on tree growth. Interestingly, the linkage between BAI and SPEI were stronger in stands with lower density (Table 1 and Fig. 4). Previous studies have pointed out that the sensitivity of Aleppo pine growth to drought is independent of stand productivity (Del Castillo et al. 2015;Del Río et al. 2014). Other studies have found geographical variation in Aleppo pine growth response to drought ) which also varies between different age classes (Martínez del Castillo et al. 2018). In uneven-aged stands, the number and characteristics of the tree samples may have profound impacts on the estimation of growth based on ringwidth measures (Nehrbass-Ahles et al. 2014;Babst et al. 2018), and probably on growth responses to drought (Martínez del Castillo et al. 2018). The analysis of even-aged plantations reduces the confounding effects of dissimilar tree level features, clarifying the impact of stand structure on climate response. Competition-mediated impacts of drought on growth have been reported in comparable Aleppo pine stands (Helluy et al. 2020). It is not surprising that reduced stand density will lead to higher growth rates and higher responsiveness to water availability, as it has been found to benefit water use efficiency (Calev et al. 2016;Tsamir et al. 2019;Manrique-Alba et al. 2020). Olivar et al. (2014) found that the removal of more than 30% of the initial basal area resulted in growth releases of Aleppo pine trees in experimental forests. Further, other factors not considered explicitly in this study may play an important role explaining tree growth. In this respect, soil attributes can modify the growth and productivity (Bravo and Montero 2001;Bravo-Oviedo and Montero 2005), although the impact of these factors might be partially considered in the site index.  It has been previously found that secondary growth measures and remotely sensed data products may correlate poorly under some circumstances (Seftigen et al. 2018;Kannenberg et al. 2019). Here we demonstrate that interannual variations in summer NDVI and secondary growth are positively related (Fig. 3) probably due to their common response to drought as found in previous studies (e.g. Xu et al. 2019). However, we found no clear convergences on how stand characteristics modulate NDVI and BAI. We found that long-term NDVI (1986-2016 period) varied from site to site as a function of stand basal area and these results point to the importance of canopy conditions on controlling primary productivity. Thus, the prevailing positive trend of NDVI in the studied sites together with its dependency on stand basal area suggest that canopy closure and densification can be the main factors driving NDVI in the long-term. Being a measure of light reflectance, NDVI will be strongly determined by the structure and composition of the vegetation. This means that stand basal area and density modify NDVI through its impact on canopy closure (Erasmi et al. 2021). All together, these results suggest inter-annual temporal couplings between BAI and NDVI reflect the sensitivity of Aleppo pine growth to drought, but that the differential sensitivity to site and stand features drives their spatial decoupling. This aligns with our hypothesis that in the complex process of tree growth, several steps are implied and thus NDVI and tree-ring width measures can provide complementary information on forest response to drought. Unfortunately, we do not have temporal direct measures of changes in canopy closure and considering the effect of changes in stand density might be useful to support our hypothesis.
We found that both, BAI and summer NDVI were reduced in response to severe droughts (see Figs. 2, 6), but these reductions were stronger in BAI as previously found in large-scale studies (Gazol et al. 2018). In terms of BAI, the 2005 drought was characterized by longer recovery periods and stronger cumulative drought impacts, hinting on the existence of lagged growth responses in the years following drought occurrence at least in this drought (i.e. legacy effects sensu Anderegg et al. 2015;see Fig. 6). In this respect, around 50% of the studied trees did not recover pre-drought BAI levels in the four years after drought occurrence. Albeit Aleppo pine is a species adapted to drought occurrence (Esteban et al. 2010;Del Río et al. 2014;Taïbi et al. 2017;Gazol et al. 2017) mortality events have been observed in plantations as well as in isolated and growth limited populations (Sarris et al. 2007;Sánchez-Salguero et al. 2010). The existence of many trees that did not recover pre-drought growth level may indicate the vulnerability of some stands to extreme droughts (De Soto et al. 2020). In terms of NDVI, results were less clear as recovery periods were shorter in 2012, but cumulative growth impacts were larger. Overall, NDVI seems to recover faster than secondary growth suggesting that canopy formation is prioritized after drought occurrence (Kannenberg et al. 2019) and potentially explaining why secondary growth is more affected by drought than NDVI (Gazol et al. 2018).
We found that short-term growth responses to drought were independent of stand variables, which contrast with results found in other semi-arid forests (Andrews et al. 2020). The lack of influence of stand variables on short-term growth responses to drought suggests that when droughts are very severe their impact on growth cannot be modulated by stand characteristics. In a future context of more frequent droughts due to climate warming, reducing tree density has been included in the toolbox of forest adaptive management practices to increase tree resistance and resilience to drought (Bottero et al. 2017;Andrews et al. 2020;Beckmann et al. 2021;D'Amato and Palik, 2021), even if the effect of this practice may not be universal (Bosela et al. 2021).

Conclusions
With the ongoing climate change, the likelihood of hotter drought occurrence, because of the concomitant occurrence of dry years and hot summers, and their impacts on forest vigour increases (Allen et al. 2015). In the Mediterranean region, warming trends are far higher than average, and dry years occur frequently, thus threatening the maintenance of natural forests and plantations (Cramer et al. 2018). Understanding which forests and under what circumstances are more vulnerable to drought is thus fundamental to mitigate the impacts of climate change. The results of this study suggest that the response of Aleppo pine growth to drought depends on stand density. It is plausible to think that management activities can thus modulate how trees within stands respond to SPEI in normal years. However, we found that stand density has no influence on the resistance of Aleppo pine to drought and its capacity to recover after it. This suggests that less dense stands have higher potential to respond to water availability that denser stands, but under severe droughts this beneficial effect might be erased. A question that emerges here is whether less dense stands may show lower mortality in response to extreme droughts. Under the ongoing scenario of more frequent severe drought events, further studies should consider this aspect.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10342-022-01526-9. data. We thank the editor and two reviewers for their constructive suggestions.
Authors contribution VR, JMO and FRP conceived of the presented idea and designed the study, VR and SCA performed the measurements, CG led the remote sensing data, AG and RAP analysed the data and led the writing, All authors discussed the results and commented on the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This work was supported by projects SEÑALES and INFRARED (VA026P17 and IR2020-1-UVA08; Junta de Castilla y León and UE-FEDER), spRING (CGL2017-87309-P; Ministerio de Ciencia e Innovación, Spanish Government) and BiOLi-DAR (PTQ-15-08067; Department of Economy, Industry and Competitiveness). Field data collection was funded by the FEADER under the provisions of the Rural Development Program of Aragón 2014-2020. AG is supported by the Ramon y Cajal Program of the Spanish MICINN under Grant RyC2020-030647-I, and by CSIC under grant PIE-20223AT003.
Data availability Should the manuscript be accepted, the datasets generated and/or analysed during the current study will be made available upon request.
Code availability Should the manuscript be accepted, the code used to analyse the data in R will be made available upon request together with the datasets.

Conflict of interest
The authors declare that they have no conflict of interest.
Ethics approval Not applicable.

Consent for publication Not applicable.
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/.