Tree Mortality Risks Under Climate Change in Europe: Assessment of Silviculture Practices and Genetic Conservation Networks

General Context: Climate change can positively or negatively affect abiotic and biotic drivers of tree mortality. Process-based models integrating these climatic effects are only seldom used at species distribution scale. Objective: The main objective of this study was to investigate the multi-causal mortality risk of five major European forest tree species across their distribution range from an ecophysiological perspective, to quantify the impact of forest management practices on this risk and to identify threats on the genetic conservation network. Methods: We used the process-based ecophysiological model CASTANEA to simulate the mortality risk of Fagus sylvatica, Quercus petraea, Pinus sylvestris, Pinus pinaster, and Picea abies under current and future climate conditions, while considering local silviculture practices. The mortality risk was assessed by a composite risk index (CRIM) integrating the risks of carbon starvation, hydraulic failure and frost damage. We took into account extreme climatic events with the CRIMmax, computed as the maximum annual value of the CRIM. Results: The physiological processes' contributions to CRIM differed among species: it was mainly driven by hydraulic failure for P. sylvestris and Q. petraea, by frost damage for P. abies, by carbon starvation for P. pinaster, and by a combination of hydraulic failure and frost damage for F. sylvatica. Under future climate, projections showed an increase of CRIM for P. pinaster but a decrease for P. abies, Q. petraea, and F. sylvatica, and little variation for P. sylvestris. Under the harshest future climatic scenario, forest management decreased the mean CRIM of P. sylvestris, increased it for P. abies and P. pinaster and had no major impact for the two broadleaved species. By the year 2100, 38–90% of the European network of gene conservation units are at extinction risk (CRIMmax=1), depending on the species. Conclusions: Using a process-based ecophysiological model allowed us to disentangle the multiple drivers of tree mortality under current and future climates. Taking into account the positive effect of increased CO2 on fertilization and water use efficiency, average mortality risk may increase or decrease in the future depending on species and sites. However, under extreme climatic events, our process-based projections are as pessimistic as those obtained using bioclimatic niche models.


INTRODUCTION
Sustainable forest management ultimately aims to attain a balance between society's increasing demands for forest products and services, and the long-term preservation of forest health and biodiversity. However, ongoing climate and global changes pose major challenges to this balance by affecting forest ecosystem functioning (Hanewinkel et al., 2012;Verkerk et al., 2020), the composition of communities and populations, and consequently the geographic distribution of species (Lindner et al., 2014). Therefore, there is an urgent need to understand the future ecological dynamics of forest ecosystems to guide sustainable management and conservation efforts.
Climate change interacts with other global change factors such as nitrogen deposition, atmospheric CO 2 concentration increase, 0 3 pollution, land use, and forest management. These factors may have antagonistic or partially offsetting effects on the physiological functioning of trees and the ecological dynamics of forests (Walther et al., 2002;Begon et al., 2006;Thuiller et al., 2006;Morin et al., 2008;Lindner et al., 2014;Pretzsch et al., 2014;Anderegg et al., 2015). For example, increasing atmospheric CO 2 concentration can increase water use efficiency on the one hand and photosynthesis on the other hand (see Walker et al., 2021 for a synthesis). Rising temperatures can also increase the length of the growing season (i.e., the period between leaf budburst and leaf fall; Davi et al., 2006b;Keenan et al., 2014) and consequently increase tree growth and forest productivity. However, the combination of rising temperatures and decreasing summer precipitation in some regions has been associated with increasing frequency and duration of droughts. This results in a decreasing tree growth and forest productivity (Zhao and Running, 2010;Reyer et al., 2014;Reyer, 2015), increasing risks of tree mortality (Allen et al., 2010) and increasing risk of disturbances (wildfire, pests; Seidl et al., 2014Seidl et al., , 2017. In addition, a single climatic driver can have both positive and negative effects on tree performance; this is exemplified by the role of increasing temperature on the advance of leaf phenology (Bigler and Bugmann, 2018), which is expected to increase vegetation length, but also to expose plants to higher risk of late frosts (Vitasse et al., 2014).
Correlative bioclimatic niche models (BNMs) are popular modeling tools used to predict the future risk of mortality and species distribution ranges (Urban, 2015). A major strength of BNMs is that the data needed to calibrate them (e.g., speciesspecific presence/absence) are available in large numbers and with increasing resolution in open-access databases (Duputié et al., 2014). However, most BNMs are not mechanistic and are usually not able to take into account the physiological response of trees to new environments (e.g., the combination of increasing CO 2 concentration and increasing temperatures). Furthermore, the spatial variation in soil properties could have complex effects on species favorability, depending on the interaction between the soil-related water capacity, the climate-related water availability, and the species-specific physiological vulnerability to water stress (E Silva et al., 2012). These well-known limitations of correlative BNMs may explain the contradictory projections sometimes generated with different BNMs, as for Abies alba along the southern edge of its distribution (Mairota et al., 2013;Tinner et al., 2013).
Ecophysiological process-based models (PBMs) offer the advantage to simulate vegetation functioning in response to explicit climate and soil variability, through their impacts on plant physiology (e.g., Cramer et al., 2001;Dufrêne et al., 2005). Regarding the prediction of future tree species distribution range, comparisons of PBMs vs. BNMs showed that PBMs were generally less pessimistic (Morin and Thuiller, 2009;Cheaib et al., 2012). Ecophysiological PBMs were initially developed to simulate carbon and water fluxes in forest ecosystems, but can also be used to investigate the environmental drivers and physiological processes triggering tree mortality under climate change. For instance, they are particularly suitable to disentangle the physiological processes contributing to mortality in face of drought (McDowell et al., 2008), and the intricate roles of hydraulic failure (the loss of hydraulic conductivity resulting from xylem embolism) and carbon starvation (the depletion of carbon reserve resulting from stomata closure to avoid hydraulic failure). These same models are also able to model the phenology of plants and the dynamics of cold resistance of the organs (Leinonen, 1996).
Forest dynamics models also have a long tradition of being used to support forest management. Recently, both ecophysiological PBMs or forest dynamics models incorporating ecophysiological processes have been applied to evaluate how management and climate change may interact to influence forest dynamics (Oddou-Muratorio et al., 2020). For instance, simulations studies supported that decreasing stand density through thinning can improve the resistance of forests to climate change by decreasing the competition for water amongst the remaining stems, and can be used to inform thinning intensity and frequencies (Guillemot et al., 2014;Cabon et al., 2018). However, the consequences of management practices on forest ecophysiological functioning remain largely unknown. Hence, it is important to investigate if and how current locally prescribed management systems may mitigate the impacts of climate change on the risk of mortality.
Assessing the future physiological functioning and ecological dynamics of forest stands requires not only accounting for the multiple effects of climate change, and for their variation across tree species distribution ranges, but also for the adaptive response of tree populations and the possible effects of forest management on this response. The adaptive potential of tree populations in face of a changing environment is usually assumed to be non negligible: besides tracking their ecological niche spatially through migration, tree populations can adapt in the short-run through individual physiological tolerance, and/or in the longrun through evolutionary response to environment-induced natural selection (Benito Garzón et al., 2011;Oddou-Muratorio and Davi, 2014). Genetic diversity also represents a key resource that can be harnessed by managers for resilience (Fady et al., 2016). Existing networks of dynamic Genetic Conservation Units (GCUs) for forest trees represent a central reservoir of possible future options for adaptation of forests to climate change (Supplementary Figure 1 for the GCUs distribution by species). In Europe, this network was set up to dynamically conserve forest genetic resources by preserving the ecological and evolutionary processes contributing to the adaptive potential of tree populations to environmental variations (Koskela et al., 2012). The European transnational network of forest tree GCUs contains 3,593 GCUs, representing 4,316 tree populations (http:// portal.eufgis.org/). BNM projections indicate that many GCUs may have a high risk of extirpation under climate change, thus weakening the conservation network as a whole (Schueler et al., 2014). These projections indicate that 33-65% of conservation units, mostly located in southern Europe, will be at the limit or outside the current bioclimatic niche of their constituting species by 2,100. The highest average increase in the risk of extirpation throughout the network can be expected for coniferous trees.
The assessment of tree mortality risk under different climate change scenarios at continental scale is needed to support sustainable forest management in Europe in three ways : (i) to identify at continental scales areas of higher risk of mortality and better understand the main drivers of mortality risk in these areas; (ii) to assess the impact of locally prescribed silviculture on tree mortality risk; (iii) and to characterize the level of threat across the in-situ gene conservation networks. Here, we focused on the risk of mortality associated with droughts and frosts, two major risks likely to increase in the future (Augspurger, 2013;IPCC, 2014;Charrier et al., 2018). We used the ecophysiological and biophysical PBM CASTANEA , to simulate the risk of hydraulic failure, carbon starvation, and frost damages of five tree species representative of the main European forest biomes: Fagus sylvatica and Quercus petraea for temperate deciduous broadleaved forests; Picea abies and Pinus sylvestris for high-latitude and high-altitude evergreen conifer forests; and Pinus pinaster, for low-latitude, evergreen conifer temperate forests (Supplementary Figure 1). Using a composite risk index of mortality (CRIM), we addressed the following issues: (1) How does the risk of mortality in face of drought and frost vary among species and across species' distribution ranges? We expected CRIM differences between species to emerge from their different sensitivity to the target climatic stress, which is modeled in CASTANEA through species-specific parameters. (2) What is the impact of forest management practices on this risk? We expected the potential mitigation effect of thinning on CRIM to vary among species and among climatic scenarios. (3) Will GCUs' risk of extinction increase under future climates, and are the projections of our PBM consistent with those based on BNMs? We expected a less pessimistic prediction of mortality risks and a lower number of GCUs with a high risk of extinction.

CASTANEA Model
CASTANEA is an ecophysiological PBM used to simulate carbon and water fluxes in forest ecosystems . Briefly, this model simulates the development of an average tree modeled by six functional compartments: canopy, branches, stem, coarse roots, fine roots, and reserves (an unlocated compartment corresponding to the Non-Structural Carbohydrates, NSC). The canopy is divided into five layers of leaves. Photosynthesis is hourly computed for each canopy layer using the model of Farquhar et al. (1980), analytically coupled to the stomatal conductance model proposed by Ball et al. (1987). The effect of temperature on photosynthesis is modeled using a response function of Rubisco-limited photosynthesis (Bernacchi et al., 2001). Maintenance respiration is modeled as proportional to the nitrogen content of the considered organs (Ryan, 1991). Growth respiration is computed from growth increment combined with a construction cost specific to the type of tissue (De Vries et al., 1974). Transpiration is also calculated hourly using the equations of Monteith (1965). The dynamics of Soil Water Content (SWC) is calculated daily using a threelayer bucket model. Soil drought drives stomata closure with a linear decrease when relative SWC is under 40% of field capacity (Sala and Tenhunen, 1996;Granier et al., 2000). In the carbon allocation sub-model Cailleret et al., 2017), the allocation coefficients of biomass between the six compartments are daily calculated, depending on the sink force and on the phenology constraints. CASTANEA was originally developed and validated at stand scale for F. sylvatica, and then for Pinus sp. and Q. petreae (Davi et al., , 2006a. The whole parameterization for the five species is detailed in Supplementary Appendix 1. In this study, we focused on five output variables simulated by CASTANEA: (1) the Net Primary Production (NPP) and ring widths to assess the ability of the model to reproduce the bioclimatic niche, (2) the percent loss of conductivity (PLC) as an indicator of risk of hydraulic failure, (3) the NSC content as an indicator of risk of carbon starvation, and (4) the frost damage index (FD).
The CASTANEA version we used allows PLC to be computed based on daily midday water potential and species vulnerability curve to embolism (Petit-Cailleux et al., 2021; Supplementary Appendix 2). To simulate budburst, we used the UniChill model (Chuine et al., 1999) in its version described in Gauzere et al. (2017). The UniChill model is a sequential two-phases model describing the cumulative effect of chilling temperatures on bud development during the endodormancy phase (first phase) and the cumulative effect of forcing temperatures during the ecodormancy phase (second phase). We simulated damages due to frosts based on the Leinonen (1996) model. Briefly, this model evaluate the daily frost damage as: Frontiers in Ecology and Evolution | www.frontiersin.org where FH is the frost hardiness (see Supplementary Material 2) and FS the frost sensitivity (fixed here to −2).
Then FD daily is summed over the year.
It is therefore sufficient for FD daily to reach 1 on a single day to cause maximum annual damage. Note that we considered that trees were able to reflush after late frosts We simulated stand mortality due to competition by using the relative density index (RDI), derived from the self-thinning rule (Reineke, 1933), as a threshold to limit stand density. RDI is classically defined as the ratio of actual stand density (N) to the maximum stand density attainable in a stand with the same average tree volume (Nmax): RDI = N/Nmax, where Nmax is defined as: with Dg the mean quadratic diameter, a the slope and b the intercept of the self thinning equation (Charru, 2012), these specific parameters are given in the Supplementary Appendix 1.
When N>Nmax, self-thinning occurs and stand density decreases to Nmax, which directly affects the Leaf Area Index (LAI) of the canopy. Indeed, the stand leaf area index (LAI stand ) is simply modeled from stand density (N), the average Leaf Mass per Area (LMA mean ), and the allometric relationship between tree size (dbh) and biomass of leaves (B leaves ).
where the biomass of leaves depends on tree diameter. Moreover, the reduction of stand density also decreases the biomass of the trunk, branches, reserves, and large roots of the average tree. All these effects (decreasing N, LAI, and biomass) were also applied where stand density is reduced by thinning (see section Management Practices below).

Climate Data and Scenarios
We considered the European area included within longitudes ranging from −11 • W to 40 • E and within latitudes ranging from 36 • N to 66 • N (continental Europe, excluding the arctic circle). CASTANEA requires the following daily climatic input variables: the minimum, mean and maximum temperatures (in • C), precipitation (mm), the wind speed (m.s −1 ), the mean relative humidity (%) and the global radiation (MJ.m −2 ). These data were derived for seven climatic scenarios (three current and four future scenarios) as detailed below (Supplementary Table 1

Current Climate
We used the Water and Global Change (WATCH)-Forcing-Data-ERA-Interim data set (WATCH in the following) to obtain current climate data at European scale (Weedon et al., 2014). This daily meteorological forcing dataset is available for the period 1979-2008 worldwide, with a spatial resolution of 0.5 • per 0.5 • . This resolution (the coarsest among the climate and soil data set) was used to divide Europe into 3,411 raster cells.

Future Climates
To take into account the uncertainties on future climatic scenarios (McSweeney et al., 2015), we used a combination of two daily regional circulation models (RCMs) developed under the EURO-CORDEX initiative (Jacob et al., 2014), with the past CO 2 concentration data and two future scenarios of representative concentration pathways (RCPs, Moss et al., 2010). All these daily meteorological model datasets are available from 1961 to 2005 (for current climatic scenario) and from 2005 to 2098 (for future climates scenarios). We selected the EUR-11.SMHI.CNRM-CERFACS-CNRM-CM5 (CM5 in the following) as warm and wet RCM (optimistic) and EUR-11.SMHI.MOHC-HadGEM2-ES (HadGEM in the following) as warm and dry RCM (pessimistic). Regarding future RCPs, the RCP4.5 scenario considers an increase of CO 2 concentration of 650 ppm with a 1.0-2.6 • C increase by 2100, and corresponds to the SRES B1 scenario (Nakićenović and Swart, 2000). The RCP8.5 scenario considers an increase of CO 2 concentration of 1,350 ppm CO 2 with a 2.6-4.8 • C increase by 2100, and corresponds to the A1F1 SRES scenario (van Vuuren et al., 2011;Harris et al., 2014;IPCC, 2014). To avoid potential biases due to differences between WATCH and EURO-CORDEX datasets, we compared the simulations run on each future EURO-CORDEX scenario to their current EURO-CORDEX reference. The two current climate datasets (CM5_hist and HadGEM_Hist) and the four future climate datasets (CM5_RCP4.5, CM5_RCP8.5, HadGEM_RCP4.5, and HadGEM_RCP8.5) were corrected for bias and downscaled using the R package "meteoland" (De Caceres et al., 2018), and considering the WATCH dataset as reference data.

Soil Data
To account for the variability of soil water capacity across Europe, we used (1) the European Soil Database to obtain data on the soil depth reached by the roots; (2) the SoilGrids 250 m database to obtain data on bulk density and clay, silt, sand and coarse fragments contents; and (3) the 3D soil hydraulic database to obtain data on soil water content at field capacity and at wilting point (Hiederer, 2013;Hengl et al., 2017;Tóth et al., 2017). All these data were aggregated from 1 × 1 km resolution to 0.5 • × 0.5 • (WGS84) resolution using the R package "raster" (Hijmans et al., 2015). Then, we extracted the mean value of each parameter at each grid point. A summary table of the values of the climate and soil variables is available in the Supplementary Table 1.

Management Practices
Four silvicultural systems (SS) were modeled on the basis of locally prescribed management practices inventoried in (Härkönen et al., 2019): • SS1, "No management": no thinning nor regeneration cuts were applied. • SS2, "Even-aged forest management with shelter-wood": intermediate cuts depend on tree height, stand basal area, and/or age. The last thinning is a shelter-wood cut, removing the overstory and coinciding with the recruitment of young seedlings (typically 10 years old). The age of the shelter-wood and intermediate cuts are specific of each species and each ecoregion. Trees that do not meet the criteria are not cut. • SS3, "Even-aged forest management with clear-cut": intermediate cuts depend on tree height, stand basal area, and/or age. After clear-cut, a new stand is planted the following year. Simulation starts when the trees reach breast height, which takes different lengths of time depending on the ecoregion and the species. • SS4, "Short rotation": there are no thinnings applied, and an integral final cut is done at an early age, followed by planting.
These silvicultural systems were only applied to the main target tree species, and were assumed to be constant over the simulated period. Thinning rules were adapted from Härkönen et al. (2019) and determined the reduction in stand density according to stand age, and eventually to stand height and/or basal area (Supplementary Table 2 Table 3).
For some countries, the share values of silvicultural systems were not available, and computed as the average of the percentage of each silvicultural systems from the closest ecoregion to the country. Supplementary Figure 4 shows the distribution of silvicultural systems across the simulated area.

Simulation Design
We first considered three past climatic scenarios (WATCH, CM5_Hist, HadGEM_Hist) and four future climatic scenario (CM5_RCP4.5, CM5_RCP8.5, HadGEM_RCP4.5, HadGEM_RCP8.5), without management. For each species and climatic scenario, we ran one CASTANEA simulation at each of the 3,411 grid points (i.e., inside and outside the current species distribution range). At each grid point, one average, even-aged stand was simulated (hence, one average tree), with fixed soil properties. For all simulations, we considered the same average starting tree, with an initial age of 8 years and a diameter at breast height of 5 cm. All species-specific parameters are available in Supplementary Appendix 1. In total, we ran 119,385 CASTANEA simulations (5 species × 7 climates × 3,411 grid points) without management (Supplementary Table 4). Secondly, we considered the management practices locally prescribed for each species. Additional scenarios with silviculture were simulated only at the grid points where at least another silvicultural systems than SS1 was identified (Supplementary Table 3). At each of these grid point, up to three additional CASTANEA simulations were run (if SS2, SS3, and SS4 were all locally applied), with the same conditions as above. In total 71,695 additional simulations with silviculture were run (Supplementary Table 4).
The entire set of simulated points under WATCH climate was used to evaluate model predictions (see next section). The simulated points within current observed species range were used to investigate the risk of mortality and its variation between climate and silviculture scenario (Supplementary Figure 1shape of species distribution is obtained from http://www. euforgen.org/species/). Note that the length of the simulated period varied among current and future scenarios (30 years for WATCH, 45 years for CM5_hist, and HadGEM_Hist, and 94 years for CM5_RCP4.5, CM5_RCP8.5, HadGEM_RCP4.5, and HadGEM_RCP8.5). Moreover, the silvicultural cycle of SS2 and SS3 most often require one hundred years to be completed. To be able to check the influence of climate change and silviculture on mortality risk and not an age effect, we compare the LAI, DBH, and stand density between WATCH, historical climate, and future climate at 2050 and 2100 years.

Model Evaluation
CASTANEA has already been evaluated at stand scale for the five species (Davi et al., 2006a;Delpierre et al., 2012). We re-evaluated the model on a series of CO 2 flux measurements at five sites. Three of these datasets were obtained from the FLUXNET2015 (Pastorello et al., 2020), database CZ-BK1 for P. abies; FR-LBr for P. pinaster; NL-Loo for P. sylvestris. The two others came from the ICOS dataset, FR-Hes for F. sylvatica and FR-Fon for Q. petraea, see Supplementary Appendix 3 for details.
To evaluate the validity of CASTANEA large-scale simulations, we used a similar approach as Cheaib et al. (2012) using "R" (R Development Core Team, 2013). We investigated the model's ability to predict species range based on thresholds values of NPP and ring widths (rw), below/above which the species is considered as absent/present. The optimal thresholds minimizing the difference between the predicted and observed distribution were selected using an iterative procedure and three statistics used in SDM evaluation: the Area Under the receiver operating Curve (AUC) and the True Skill Statistic (TSS), two classical statistics (Fourcade et al., 2018); and also the Boyces index (BI), another statistic recently proposed and non correlated with classical metrics (Konowalik and Nosol, 2021).
First, AUC represents the rate of true positives (fraction of presences that are actually detected) vs. the rate of false positives (fraction of absences that are incorrectly detected as presence). It is independent from the threshold value used to convert probabilities of presence per pixel into presence-absence data (Elith and Graham, 2009). An AUC-value of 1 means a perfect fit between the predicted and observed distribution; an AUC-value of 0.5 corresponds to random prediction of a species presence and an AUC-value value of 0 means that the prediction is always wrong. To calculate the AUC we used the "PresenceAbsence" package (Freeman, 2012). Second, the TSS compares the number of correct predictions, minus those attributable to random guessing, to that of a hypothetical set of perfect predictions (sensitivity plus specificity −1) and varies between −1 (poor fit) and 1 (perfect fit). TSS normalizes the overall presence-absence by the presence-absence that might have occurred by chance alone and is not influenced by the distribution size. To calculate TSS, we converted the net primary production and ring width simulated into presence-absence data by using a threshold value maximizing sensitivity and specificity as recommended by Liu et al. (2013). Finally, the BI only requires presences and measures how much model predictions differ from random distribution of the observed presences across the prediction gradients (Boyce et al., 2002). It varies between −1 and +1. Positive values indicate that the model's predictions are consistent with the distribution of presences, values close to zero mean that the model is not different from a random model, negative values indicate counter predictions, i.e., poor ability to predict areas where presences are more frequent (Hirzel et al., 2006). We used the "ecospat" package to calculate the BI (Broennimann et al., 2021).
Finally, we also evaluated the effects of CO 2 enrichment under future climate of the main output variables simulated by CASTANEA. To quantify and standardize CO 2 effects across variables and varying ranges of [CO 2 ] (Walker et al., 2021) proposed a relativized β factor: where Y a and Y e are the values of any response variable at lower CO 2 concentration (CO 2 ,a) and higher CO 2 concentration (CO 2 ,e), respectively. We computed the beta growth primary production (β GPP ), β WUE , and β BG simulated by CASTANEA for F. sylvatica, using simulations between 1975 and 2005 under the HadGEM_Hist scenario (ambient CO 2 ) on the one hand, and between 2070 and 2098 under the HadGEM_RCP8.5 scenario (elevated CO 2 ) on the other hand. Then, we compared these simulated values to the observations compiled in Walker et al. (2021).

Computation of the Risk of Mortality
Note that we did not directly simulate mortality due to drought and frost damage with CASTANEA because the thresholds in PLC, NSC, and FD triggering mortality are unknown. Instead, to evaluate the risk of mortality, and compare it between scenarios and species, we computed the relative values of three CASTANEA output variables over all the simulated period, as in Petit-Cailleux et al. (2021).
First, the relative percent loss of conductivity (rPLC) was computed as a proxy of the risk of mortality due to hydraulic failure: where PLC mean is the mean of yearly PLC-values over the simulated period, and PLC species threshold is the species-specific PLC-value above which mortality occurs (88% for deciduous and 50% for coniferous species; Brodribb, 2009;Urli et al., 2013). Hence, rPLC varied between 0 and 1, and increased with the risk of hydraulic failure. Because rPLC is computed from the mean PLC value over the simulated period, it is expected to buffer the effect of extreme climate events such as intense drought. To overcome this limitation, we also computed rPLC max by considering instead the maximum PLC-value over the simulated period (i.e., replacing PLC mean by PLC max in Equation 3).
Second, the relative non-structural carbohydrate content (rNSC) was computed as a proxy of the risk of mortality due to carbon starvation: if NSC mean > NSC species threshold (7) where NSC frac n is the average fraction of NSC biomass over the biomass of all the other compartments at a given year n, and NSC species threshold the threshold in NSC-value above which mortality through carbon starvation is unlikely. NSC species threshold was arbitrarily set at the CASTANEA initial NSC concentration value (g C .g C −1 ) that corresponds to the average measured value for each considered species (see references in Supplementary Appendix 1). Hence rNSC varied from 0 to 1, increasing values indicating an increasing risk of carbon starvation.
We also computed rNSC max considering instead the minimum NSC-value over the simulated period (i.e., replacing NSC mean by NSC min in Equation 4). Taking the average measured NSC value as a threshold value below which starvation risk occurs may seem too severe, but in reality, when carbon starvation occurs, the NSC drop quickly and we verified that the choice of the threshold had little impact.
Thirdly, the level of frost damage (rFD) was computed as a proxy of the risk of mortality due to frosts. The FD yearly index described in equation (2) is already relative and varies from 0 to 1, increasing values indicating an increasing risk of frost damage. rFD was thus computed as the mean of FD yearly values over the simulation period, and rFD max as the maximal FD value over the simulation period.
Note that for scenarios with silviculture, there were possibly up to four simulations, and hence four values of PLC, NSC, or FD at a given grid point. In these cases, we computed a ratio taking into account the share of each silvicultural systems j at each grid point i: where share ij is the share of each silvicultural system j at grid point i, and N is the number of different silvicultural systems at point i. The different PLC, NSC, or FD values for each silvicultural system j were weighed by Ratio ij . Finally, we computed a composite risk index of mortality (CRIM) combining rPLC, rFD, and rNSC with an identical weight. Annual CRIM values were computed as: The mean CRIM value was then averaged across the simulated period. CRIM can vary between 0 (minimal risk of mortality) and 1 (maximal risk of mortality) and CRIM n =1 as soon as either rFD max , rPLC or rNSC = 1. A similar CRIM max value was calculated by combining rPLC max , CRIM max , and rFD max . Both the Net Primary Production (NPP) and ring width (rw) simulated by CASTANEA were used to predict species distribution, based on a threshold value above/below which the species was predicted as present/absent. To select this threshold value and evaluate the quality of the prediction we used three statistics (AUC, TSS, and BI, see text for details).

Model Evaluation
We obtained a high, positive correlation of GPP, the net ecosystem exchange and ecosystem respiration between measures and simulations for broadleaved species and P. sylvestris (r ≥ 0.74; Supplementary Appendix 3.2). For P. abies and P. pinaster, the correlation between measures and observations was high for GPP and the ecosystem respiration, but low for net ecosystem exchange. We determined a species-specific threshold for ring widths and NPP-values maximizing the agreement between observed and predicted distributions for the TSS calculation. The current distribution ranges of broadleaved species were predicted with high TSS-values, high AUC-values and high BI-values; for F. sylvatica (TSS ≥ 0.62, AUC ≥ 0.86, and BI ≥ 0.94); and for Q. petraea (TSS ≥ 0.69, AUC ≥ 0.92, and BI = 1; Table 1). By contrast, the distribution ranges of coniferous species were poorly predicted: for P. abies TSS ≤ 0.15, AUC ≤ 0.56, and BI ≤ −0.03; and for P. sylvestris TSS ≤ 0.07, AUC ≤ 0.61 and with a negative correlation according to BI ≤ −0.31. For P. pinaster, the evaluation is contrasted between indices, with a poor TSS ≤ 0.32 and AUC ≤ 0.64, but a rather good prediction according to BI ≥ 0.64.

Variations of the Risk of Mortality Over Species Realized Range Under Current Climate
Under the WATCH climatic scenario, the mean CRIM simulated from 1979 to 2008 was 0.13 for F. sylvatica, 0.25 for Q. petraea, 0.25 for P. pinaster, 0.28 for P. abies, and 0.34 for P. sylvestris, within their respective realized niche (Figure 1,   Supplementary Table 5). For broadleaved species, the highest CRIM-values occurred in the eastern Mediterranean region, in mountainous areas and in the eastern part of their distribution. For P. abies and P. sylvestris, the highest CRIM-values occurred in the northern and eastern parts of their distribution. For P. pinaster, the highest CRIM-values occurred in the southern part of its distribution. When considering CRIM max to better account for the effect of extreme climatic events on the risk of mortality, mean CRIM max values were much higher for all species: 0.69 for P. pinaster, 0.48 for F. sylvatica, 0.83 for P. abies, 0.72 for Q. petraea, and 0.84 for P. sylvestris (Figure 2,  Supplementary Table 5). The spatial variations of CRIM maxvalues were similar to that of CRIM-values, except that high CRIM max values were also observed at the southern part of P. abies and P. sylvestris distributions (Figure 2).
The processes contributing to the risk index of mortality varied among species and spatially (Figures 1, 2,  Supplementary Table 5). For Q. petraea and P. sylvestris, CRIM was mostly driven by the risk of hydraulic failure (rPLC = 0.18 for both species), and CRIM max by both the risk of hydraulic failure (rPLC max = 0.6 and 0.63, respectively) and frost damage (rFD max = 0.52 and 0.57, respectively). The highest rPLC-values occurred in the Mediterranean region and in Eastern Europe for these two species (Figures 1, 2). For P. abies, the risk of frost damage drove both the CRIM-value (rFD = 0.19) and the CRIM max -value (rFD max = 0.76). The highest rFD-values occurred in Northern and Eastern Europe for this species. For P. pinaster, the risk of carbon starvation drove the CRIM-value (rNSC = 0.13), while the three risks contributed equally to CRIM max . For this species, the highest rNSC-values occurred along the coastline. Finally for F. sylvatica, CRIM was driven both the risk of hydraulic failure (rPLC = 0.07) and late frost (rFD = 0.05), while the three risks contributed equally to CRIM max (CRIM max = 0.21, rPLC max = 0.23, and rFD max = 0.30). For a given risk, the highest values did not occur in the same place for different species. For instance, for broadleaved species, the highest rFD-values mainly occurred at high altitudes (Meseta, Alps, Carpathian mountains), while for the coniferous species, the highest high rFD-values occurred in northern Europe. For a given species, the contributions of the physiological processes to the risk of mortality varied across its range. For instance considering CRIM max and Q. petraea, in the eastern continental part of the range all three risks were high, in the Alps the risk of frost was predominant, in Northern Europe (Finland) the risks of hydraulic failure and frost damage were the highest, while in the Southern part of the range (Iberian peninsula) the risk of carbon starvation was predominant.

Evolution of the Risk of Mortality Between Current and Future Climates
The variation of the risk of mortality (estimated by CRIM and CRIM max ) between current (CM5hist and HadGEM_hist) and future (CM5_RCP4.5, CM5_RCP8.5,HadGEM_RCP4.5, HadGEM_RCP8.5) climatic scenarios differed between species (Figures 3, 4, Supplementary Figure 5). In average over all pairwise comparisons, the risk of mortality increased for P. FIGURE 1 | Spatial distribution of mean values of mortality risk indices under WATCH current climate, across species' distribution range. The composite risk index of mortality (CRIM), given in the first column, is the combination of the three indices displayed on the other columns: the frost damage index (rFD) as a proxy of the risk of mortality due to frost; the relative non-structural carbohydrate content (rNSC) as a proxy of the risk of carbon starvation; the relative percent loss of conductivity (rPLC) as a proxy of the risk of hydraulic failure. All indices vary from 0 (purple, low risk) to 1 (yellow, high risk).
The maximum value of CRIM max = 1 was reached in many locations particularly for all conifers (Figures 2, 4). The risk of mortality estimated with both CRIM and CRIM max tended to be lower for the "optimistic" CM5 RCM compared to the "pessimistic" HadGEM RCM for the coniferous species. For all species except P. abies, the CRIM-values were similar for the FIGURE 2 | Spatial distribution of maximum values of mortality risk indices under WATCH current climate, across species' distribution range. The risk index of mortality (CRIM max ), given in the first column, is the combination of the three indices displayed on the other columns: the frost damage index (rFD max ) as a proxy of the risk of mortality due to frost; the relative non-structural carbohydrate content (rNSC max ) as a proxy of the risk of carbon starvation; the relative percent loss of conductivity (rPLC max ) as a proxy of the risk of hydraulic failure. All indices vary from 0 (purple, low risk) to 1 (yellow, high risk). RCP 4.5 and RCP 8.5 scenarios whereas, for all species, CRIM max tended to be higher for the RCP 8.5 scenarios (Figure 4).
To verify whether the observed changes in risks between present and future climatic scenarios could be due to the shorter simulation period for the present climatic scenario (45 years) than for the future (94 years), we also analyzed the mortality risk under future climate at the 2050 date (i.e., over a 45-year period): for each climate change scenario, the risks predicted at this intermediate date did not differ in any major way from those based on the entire simulation period (Supplementary Table 5, Supplementary Figures 6-9). Mean DBH slightly increased and tree density slightly decreased in 2050, and these changes were much more pronounced in 2100 for broadleaved species only (Supplementary Figures 10, 11).

Effect of Silviculture on the Risk of Mortality
In average, silviculture had contrasting effects on the risk of mortality under current climate depending on the species (Figure 5, Supplementary Table 5). Under HadGEM current climate, CRIM and CRIM max decreased by 0.02 and 0.01, respectively for P. sylvestris between scenario SS1 (no sylviculture) and scenarios mixing SS1-SS4 (with silviculture). This decrease was driven by a decreased risk of carbon starvation ( rNSC = −0.12; rNSC max = −0.13), despite an increased risk of hydraulic failure. However, CRIM and CRIM max increased by 0.05 and 0.04, respectively for P. abies. Silviculture had no effect on the risk of mortality for P. pinaster, F. sylvatica, and Q. petraea (Figure 5, Supplementary Table 5). Under current climatic conditions, the local impacts of management practices on the different risks were globally uniform and low all over species' range, except for P. sylvestris in which case silviculture tends to reduce the risk of carbon starvation in the northeastern part of the range, i.e., precisely where this risk is higher (Supplementary Figure 12).
Applying today's silviculture rules under future climatic scenario resulted in increasing CRIM-values as compared to the FIGURE 4 | Variation of the mean values of mortality risk indices (CRIM max , rFD max , rNSC max , and rPLC max ) under two current and four future climatic scenarios. Each future climatic scenario (in blue for the CM5 RCM and in orange for the HadGEM RCM) can be compared to their historical reference (in light and dark gray for the CM5 and HadGEM RCMs, respectively). Simulations were run without management (SS1). See legend of Figure 3 for whisker boxes.
baseline scenario without silviculture (hist) for P. abies (mean CRIM = 0.04) and P. pinaster (mean CRIM = 0.03). This increase in CRIM was mostly driven by increasing risk of hydraulic failure (mean rPLC = 0.03 and 0.01 for P. abies and P.pinaster, respectively) and carbon starvation (mean rNSC = 0.01 and 0.04 for P. abies and P. pinaster, respectively). By contrast, silviculture was always associated to decreasing CRIM and rNSC-values for P. sylvestris (mean CRIM = 0.10 and mean rNSC = −0.22). However, CRIM max -values increased for all species, including P. sylvestris, between scenarios without vs. with silviculture under future climates (mean CRIM max = 0.05 over the five species). This increase in CRIM max was mostly driven by increasing risk of carbon starvation (mean CRIM max = 0.15 over the five species) (Figure 6, Supplementary Table 5).
FIGURE 5 | Variation of the mean values of mortality risk indices (CRIM, rFD max , rNSC, and rPLC) under two current (CM5-hist in light gray and HadGEM-hist in dark gray) and two future (CM5-8.5 in blue and HadGEM-8.5 in orange) climatic scenarios. Each scenario with locally prescribed silviculture (i.e., four possible silvicultural systems) is represented with gray whisker box contouring, and can be compared to its reference without management (i.e., SS1), with black contouring. See legend of Figure 3 for whisker boxes.
These changes in CRIM max were associated to changes in tree density and LAI. In average, silviculture (SS1-SS4) was associated to lower LAI and density per hectare in P.abies and P.sylvestris. Silviculture was associated to lower LAI only for F. sylvatica and to lower density per hectare only for P. pinaster, while it has no effect on LAI and density per hectare of Q. petraea (Supplementary Figures 11, 13, 14).

The Risk of GCU's Extirpation
We used the value of CRIM max = 1, i.e., predicted occurrence of mortality of the simulated stand during a period of time, as a criterion of threat of extinction. Under the most severe climate change scenario (HadGEM_rcp8.5), from 38% (F. sylvatica) to 90% (P. pinaster) of the GCUs are at threat of extinction in 2100 assuming that they are not managed, and these values are slightly FIGURE 6 | Variation of the mean values of mortality risk indices (CRIM max , rFD max , CRIM max , and rPLC max ) under two current (CM5-hist in light gray and HadGEM-hist in dark gray) and two future (CM5-8.5 in blue and HadGEM-8.5 in orange) climatic scenarios. Each scenario with locally prescribed silviculture (i.e., four possible silvicultural systems) is represented with gray whisker box contouring, and can be compared to its reference without management (i.e., SS1), with black contouring. See legend of Figure 3 for whisker boxes.
higher if we assume they are managed following the current local silvicultural practices ( Table 2), from 40% (F. sylvatica) to 98% (P. pinaster) in this case. The evolution of the index of mortality risk (CRIM max ) between current and 2100, under the harshest climate change scenario (HadGEM_rcp8.5), is highly variable among species: for F. sylvatica the risk will increase in 19% and decrease in 60% of the GCUs, whereas for P. pinaster the risk will increase in 48% and decrease in 2% of the GCUs ( Table 2).

DISCUSSION
In this study, we used the process-based model CASTANEA to simulate spatial and temporal variation in mortality risks due  to hydraulic failure, carbon starvation and frost damage. These risks were assessed in five major European tree species across their distribution range, as a function of variations in climate, soil properties and management practices. We focused on the risk of mortality associated with droughts and frosts (winter, spring, and fall frosts with varying frost resistance across season), two major risks likely to increase in the future (Augspurger, 2013;IPCC, 2014;Charrier et al., 2018).

The Impacts of Drought and Frost Vary Between Species and Across Distribution Range
Considering the current climate, our results show that the physiological processes driving the risk of mortality differ among species. Carbon starvation drives the mean mortality risk (CRIM) for P. pinaster; frost damage drives the CRIM for P. abies; hydraulic failure drives the CRIM for Q. petraea and P. sylvestris; and a combination of hydraulic failure and frost damage drives the CRIM for F. sylvatica. When considering extreme climatic events with CRIM max , the risk is higher and the contributions of the three focal physiological processes are more balanced. The rank of simulated species-specific risks of hydraulic failure or frost damage were not necessarily those expected from physiological or phenological observations. Indeed, deciduous broadleaves have a narrower safety margin for hydraulic failure as compared to evergreen conifers (Choat et al., 2012;Martin-StPaul et al., 2017). Late-leafing species such as P. abies have a higher safety margin for frost damage as compared to early-leafing species (Bigler and Bugmann, 2018). The low contribution of carbon starvation to the mean risk index is consistent with current ecophysiological knowledge (Martínez-Vilalta et al., 2016). However, the fact that simulations are initiated with small trees (diameter 5 cm) probably contribute to explain this result. Simulations run with bigger and older trees on average may show higher carbon starvation risk. We also find that the physiological processes driving the risk of mortality vary across species distribution ranges. The highest risk of hydraulic failure is found under Mediterranean and continental climates in Eastern Europe, associated to long and severe droughts, which is consistent with other studies based on PBMs or cSDMs (Gárate-Escamilla et al., 2019). As expected, the places the most vulnerable to frosts are located at high latitude and altitude. Secondly, our results highlight that the evolution of mortality risk under future climate by 2100 differs between species. CRIM markedly decreases for broadleaved species and P. abies, and slightly decreases for P. sylvestris. This decrease is driven by decreasing risks of hydraulic failure and frost damage, which compensates for the increasing risk of carbon starvation for coniferous species. Broadleaved species remain at a low risk of carbon starvation under future climate, which can be due to their higher stomatal sensitivity to CO 2 . Indeed, Klein and Ramon (2019) showed that conifers will need as much water as they currently do to complete their photosynthesis under future climate, while angiosperms will need less water for the same amount of CO 2 to complete their photosynthetic cycles. In contrast, for P. pinaster, the CRIM increases due to increasing risk of frost damage and carbon starvation. For the broadleaved species and P. abies, our model predicts a marked reduction of CRIM between current climatic conditions and future RCP 4.5 scenario, but not so much changes between RCP 4.5 and RCP 8.5 scenarios. This singular pattern suggests a non linear beneficial effect of CO 2 which compensates for the other harmful effects of climate change.
The evolution of mortality risk under future climate by 2100 also differs through space. Overall, for coniferous species, the risk of carbon starvation increases everywhere except in areas where it was already high and in mountainous areas. The risk of frost damage decreases or remains stable in most of Europe, which may be due to less severe winters; however it increases for P. pinaster in the Iberian Peninsula, likely due to earlier budburst and late frost in this area (Vitasse et al., 2014). The risk of hydraulic failure tends to decrease slightly throughout Europe between the current and future climate, except in a few patches mainly due to soil characteristics.

The Impact of Current Management Practices on the Risk of Mortality Under Current and Future Climate Varies Between Species
In our simulations, considering the current management practices has an ambivalent effect on the risk of mortality. Under current climate conditions, local forest management practices decrease in average the mean risk of mortality (CRIM) for P. sylvestris, slightly increase it for P. abies and has no impact on the other species. Under future climate conditions, the same local forest management practices still decrease the CRIM for P. sylvestris, slightly increase it for P. abies and P. pinaster and has no impact for the other species. Moreover, for all species, silviculture increases the risk of mortality due to carbon starvation when considering extreme drought events (rNSC max ).
Our results are thus at odds with previous studies showing that forest thinning could mitigate drought impacts (Elkin et al., 2015;Sohn et al., 2016). The simulated management practices reduce the average density of forest stands in coniferous species. This leads to a decrease of the average tree biomass and thereby its respiration. But the leaf area index (LAI) and thus the photosynthesis are also reduced by thinning. Hence, the effect of management practices on the risk of carbon starvation simulated by CASTANEA depends on the respective effects of forest thinning on respiration on the one hand, and photosynthesis on the other hand. For P. sylvestris this effect is always stronger on respiration and therefore management practices decreases the carbon starvation risk under current and future climates. But for P. abies and P. pinaster, management practices have the opposite effect in the future and increase the risk of carbon starvation.
The leaf area index decrease should also mitigate the risk of hydraulic failure by decreasing evapotranspiration (in average by 6% in our simulations, see Supplementary Table 5). That is why the increase in rPLC with silviculture for some species is quite unexpected. In these cases, the decrease in leaf area index is greater than the decrease in transpiration (data not shown), because leaf area index is quite high. The leaf midday water potential that determines the risk of cavitation is not calculated from the total tree transpiration, but from the transpiration per m 2 of leaf. Physiologically, the higher the flux at leaf level, the greater the risk of cavitation. When the canopy is dense there are more leaves and therefore the flux is distributed over more surface, which explains why in these cases, rPLC is finally lower without silviculture.
This study is among the first attempts to simulate the largescale effects of management practices with an ecophysiological PBM (see also Härkönen et al., 2019). Our simulations suggest that the impact of silviculture on the risk of mortality might not be so obvious and depend on species, the sites and the climatic scenarios considered. In particular for broadleaved species, the impacts of silviculture on the risks of mortality are uniformly low over species' distribution range. To improve these predictions, future simulations should take into account the variations in age and size class distribution at initialization.

The Threat of GCU Extirpation Under Climate Change Differs Among Species
Under the harshest climate change scenario, our projections of mortality risk in the five in situ conservation networks predict that 38-90% of the GCUs are at threat of extinction in 2100, depending on the species. These PBM-based projections are more pessismistic than earlier BNM-based projections predicting that 33-65% of the GCUs will be at the limit or outside their current climatic niche (Schueler et al., 2014). This unexpected contrast between PBM and BNM predictions has two main explanations. First, in contrast to Schueler et al. (2014) assumption that 95% of the species' current climatic niche is considered to be riskless, our CRIM max predictions indicate that most of the CGUs are already in a high mortality risk zone under current climate. This is consistent with the policy of choosing some GCUs located in areas that are already experiencing climate-related selection pressures. There are few GCUs (from 5 to 31%) shifting from low CRIM max (<0.75) in the current climate to high CRIM max in the future climate (>0.75); by contrast, 12-51% of CGUs shift from high CRIM max to low CRIM max . Secondly, and contrasting with the line of previous studies (Morin and Thuiller, 2009;Cheaib et al., 2012), we used the maximal and not the mean risk index of mortality to predict the threat of GCU extirpation. We considered spatiotemporal patterns of mean CRIM variation useful to disentangle how hydraulic failure, carbon starvation and frost damage respectively contribute to the risk of mortality; however, we argue that realistic predictions of the extirpation threat should rather account for extreme climatic events triggering maximal physiological risks of mortality. Doing so,we may overestimate this threat. In particular, because of the coarse spatial resolution of our simulation grid (0.5 × 0.5 • ), the local climates at actual CGUs location may differ from the average climate at the grid point.
Our projections also show that traditional management systems would not help to reduce the mortality risk in those areas, and that risk reduction would require other specifically dedicated management practices. Our PBM approach lies on the pessimistic side of the projections because it does not account neither for the within-stand diversity of the response to climate, which is a factor of resilience of the populations facing disturbance since not all trees are equally threatened at the same time, nor for the evolutionary potential of each forest stand, which is a factor of adaptation from one generation of trees to the next. However, in the harshest situations, silvicultural interventions and natural evolution might not be sufficient to conserve the local population in situ, and ex situ conservation strategies, either dynamic or static, should be considered.
The EUFORGEN program has developed an operational decision support tool for the management of GCU based on demographic, genetic and disturbance indicators that can be measured by different verifiers, and recommended actions depending on the levels of these indicators (Rudow et al., 2020). A combined index of mortality risk such as CRIM max could be used as a verifier to inform the risk of population loss.

Future Direction for PBM Improvements
The growth and NPP simulated with CASTANEA across the whole Europe (i.e., potential niche) under current climate were consistent with the realized distribution range of broadleaved species, but less so for coniferous species. Three hypotheses might explain why NPP variations may less correctly predict the realized distribution of coniferous species: (1) their distribution is rather determined by biotic interactions and management practices rather than by climate and soil conditions; (2) the coarse climatic grid does not correctly represent local climate conditions in mountainous areas where P. sylvestris and P. abies are located; (3) the realized distribution is not limited by survival at adult stage, but by reproduction or survival at seedling stages (Niinemets and Valladares, 2006). In addition, for P. abies, in the comparison with the CO 2 flux measurements, we may underestimate the GPP if the average value of Vcmax used for the simulations and obtained from the literature underestimate photosynthesis (Walker et al., 2014), which can lead to too low values especially at the northern and eastern boundaries. The large-scale validation of CASTANEA needs to be improved for some species to better predict mortality. Two approaches, out of the scope of the current study, can be used. The first approach would be to simulate mortality at stand scale and compare it with mortality observations at European scale, as done in Petit-Cailleux et al. (2021) at local scale. A second one would be to compare NPP simulated with those measured from satellite data. These two approaches would additionally require to account for the variations in age and size classes for forest stand across the simulated area.
CASTANEA simulated a strong positive effect of increased CO 2 on Gross Primary Production (GPP), Water Use Efficiency (WUE) and biomass growth (BG). The simulated β GPP is close to average value obtained across FACE experiments (β GPP = 0.73, 11 sites, 45 species). The simulated β WUE is a bit lower than the observed increase in ecosystem-scale plant WUE from Duke and ORNL FACE experiments (0.76 and 1.1, respectively). Finally, the simulated β BG is higher than the estimated β BG (between 0.49 and 0.56) from FACE experiments (Walker et al., 2021). However, CASTANEA does not take into account the possible acclimation of photosynthesis to rising CO 2 concentration and temperatures neither that of respiration to rising temperatures and water stress. Yet, in C3 plants, the maximum carboxylation rate and the maximum rate of electron transport were significantly reduced at elevated CO 2 (Ainsworth and Rogers, 2007). Most species can also shift their thermal optimum of photosynthesis upward in response to warming (Crous, 2019). Similarly, temperature acclimation of respiration involves a change in respiratory capacity and/or the sensitivity of respiration rate to temperature (Atkin et al., 2008). These acclimation processes would need to be accounted for in future simulation studies.
We also only partially accounted for interspecific variability in thresholds to mortality (e.g., the PLC value leading to mortality differed among broadleaved and conifer trees), and not at all for intra-specific variability in vulnerability. Yet, local adaptation is widespread in tree populations throughout their distribution ranges (Benito Garzón et al., 2011;Alberto et al., 2013). For example, Anderegg et al. (2016) demonstrated a genetic differentiation of xylem resistance traits in several species. For F. sylvatica, Kreyling et al. (2014) found genetic differentiation of frost resistance parameters, particularly in marginal populations. These patterns of intra-specific variability and differentiation of adaptive traits involved in response to drought and frost would be interesting to consider in future simulations.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories at: doi: 10.5281/zenodo.5160130.

AUTHOR CONTRIBUTIONS
SO-M, HD, FL, and BF initiated the study. SO-M, ML, and PV reviewed the silvicultural management across Europe. CP-C, HD, and SO-M upgraded and ran the PBM. CP-C, SO-M, and FL created the post processing method. CP-C drafted the manuscript. All authors contributed to its improvement and approved its final version.

ACKNOWLEDGMENTS
This work was based on the Capsis platform, http://www. inra.fr/capsis, and we thank François de Coligny for his help. We thank Marta Benito-Garzon, François Pimont and Valentin Journé for useful discussion that helped improving the manuscript. The data set FLUXNET2015 allowed us to validate the model. Parametrization of the model have benefit from the request 6015 (Modeling vulnerability of European trees under climate change) to the TRY database from 2019-04-02.