Shifts in Forest Species Composition and Abundance under Climate Change Scenarios in Southern Carpathian Romanian Temperate Forests

: The structure and functioning of temperate forests are shifting due to changes in climate. Foreseeing the trajectory of such changes is critical to implementing adequate management practices and deﬁning long-term strategies. This study investigated future shifts in temperate forest species composition and abundance expected to occur due to climate change. It also identiﬁed the ecological mechanisms underpinning such changes. Using an altitudinal gradient in the Romanian Carpathian temperate forests encompassing several vegetation types, we explored forest change using the Landis-II landscape model coupled with the PnET ecophysiological process model. We speciﬁcally assessed the change in biomass, forest production, species composition and natural disturbance impacts under three climate change scenarios, namely, RCP 2.6, 4.5 and 8.5. The results show that, over the short term (15 years), biomass across all forest types in the altitudinal gradient will increase, and species composition will remain unaltered. In contrast, over the medium and long terms (after 2040), changes in species composition will accelerate, with some species spreading (e.g., Abies alba Mill.) and others declining (e.g., Fagus sylvatica L.), particularly under the most extreme climate change scenario. Some forest types (e.g., Picea abies (L.) karst forests) in the Southern Carpathians will notably increase their standing biomass due to climate change, compared to other types, such as Quercus forests. Our ﬁndings suggest that climate change will alter the forest composition and species abundance, with some forests being particularly vulnerable to climate change, e.g., F. sylvatica forests. As far as productivity and forest composition changes are concerned, management practices should accommodate the new conditions in order to mitigate climate change impacts.


Introduction
Human-induced climate change is one of the major processes affecting the global environment nowadays [1]. However, the impacts are so complex and diverse that the net effect of climate change on forest systems is still uncertain. For instance, while the increase in atmospheric CO 2 concentrations often leads to greater productivity [2,3], the aridity caused by warming [4,5] and the intensification of disturbances regimes [6] usually have a negative impact on forests' structure and productivity. Recent studies have suggested that despite the positive climate change-associated effects [7,8], the negative ones often prevail [9,10].

Study Area
The study area is located in the Fagaraş Mountains, in the Romanian Southern Carpathians (Figure 1). It occupies around 6400 km 2 and comprises a large geographical gradient that extends 110 km north to south and 73 km east to west, covering the main forested lands in the Argeş, Sibiu and Braşov counties. The altitude ranges from 185 to 2544 m.a.s.l, with plain landscapes (c. 12% average slope, elevation < 500 m) on the southern sector and steep and rough topography (c. 32% average slope, elevation > 1000 m) on the high-altitude northern sector. The climate varies along the altitudinal gradient [23]: the annual mean temperature in the highlands is lower than 4 • C, with a winter mean temperature below −5 • C and a summer temperature around 14 • C. The precipitation approaches 1000 mm, with summer rainfall of around 315 mm. The climate in the plains is warmer and drier, with a 10.5 • C annual mean temperature, mean winter temperatures below 0 • C and summer temperatures over 20 • C. The annual precipitation here is about 700 mm, with summer rainfall of around 260 mm. In general, September and October are the driest months, while May and June are the wettest ones. Soil types in the area include protosols, spodozols, cambisols, argiluvisols and hydromorphic soils with pseudogleic properties [34,35].

Study Area
The study area is located in the Fagaraş Mountains, in the Romanian Southern Carpathians ( Figure 1). It occupies around 6400 km 2 and comprises a large geographical gradient that extends 110 km north to south and 73 km east to west, covering the main forested lands in the Argeş, Sibiu and Braşov counties. The altitude ranges from 185 to 2544 m.a.s.l, with plain landscapes (c. 12% average slope, elevation < 500 m) on the southern sector and steep and rough topography (c. 32% average slope, elevation > 1000 m) on the high-altitude northern sector. The climate varies along the altitudinal gradient [23]: the annual mean temperature in the highlands is lower than 4 °C, with a winter mean temperature below −5 °C and a summer temperature around 14 °C. The precipitation approaches 1000 mm, with summer rainfall of around 315 mm. The climate in the plains is warmer and drier, with a 10.5 °C annual mean temperature, mean winter temperatures below 0 °C and summer temperatures over 20 °C. The annual precipitation here is about 700 mm, with summer rainfall of around 260 mm. In general, September and October are the driest months, while May and June are the wettest ones. Soil types in the area include protosols, spodozols, cambisols, argiluvisols and hydromorphic soils with pseudogleic properties [34,35]. Temperate forests, whose composition changes across the altitudinal belts, cover 52% of the total study area. There are five general forest types: (a) Picea abies (L.) karst forests, (b) mixed Fagus sylvatica L.-conifer forests, (c) F. sylvatica forests, (d) mixed broadleaved forests and (e) Quercus forests. Picea forests (Norway spruce) occupy the subalpine and montane superior belts. Mixed Fagus-conifer forests, mainly Picea and Abies alba Mill., dominate the intermediate montane belt, with the beech abundance increasing as the altitude decreases. Mixed Fagus-conifer forests are substituted by pure Fagus forests first and mixed broadleaved forest (F. sylvatica and Carpinus betulus L.) later, both of them in the inferior montane belt [24,36]. Oaks gradually appear in mixed broadleaved forests in the colline belt and, in low-altitude areas, reach the point where Quercus forests dominate the landscape. Within Quercus forests, there is an altitudinal transition from Quercus petraea (Matt.) Liebl. to Quercus robur L., Quercus cerris L. and Quercus frainetto Ten. in the most southern thermophilus areas [36]. Temperate forests, whose composition changes across the altitudinal belts, cover 52% of the total study area. There are five general forest types: (a) Picea abies (L.) karst forests, (b) mixed Fagus sylvatica L.-conifer forests, (c) F. sylvatica forests, (d) mixed broadleaved forests and (e) Quercus forests. Picea forests (Norway spruce) occupy the subalpine and montane superior belts. Mixed Fagus-conifer forests, mainly Picea and Abies alba Mill., dominate the intermediate montane belt, with the beech abundance increasing as the altitude decreases. Mixed Fagus-conifer forests are substituted by pure Fagus forests first and mixed broadleaved forest (F. sylvatica and Carpinus betulus L.) later, both of them in the inferior montane belt [24,36]. Oaks gradually appear in mixed broadleaved forests in the colline belt and, in low-altitude areas, reach the point where Quercus forests dominate the landscape. Within Quercus forests, there is an altitudinal transition from Quercus petraea (Matt.) Liebl. to Quercus robur L., Quercus cerris L. and Quercus frainetto Ten. in the most southern thermophilus areas [36]. Some of the forested lands are actively managed for timber production, with F. sylvatica providing 34% of the overall standing wood volume, P. abies 28%, Quercus spp. 17%, A. alba 5% and other species the remaining 16% [37][38][39][40][41][42][43][44][45][46]. Depending on forest species and conditions, different management systems are implemented including tree selection, shelterwood and clearcutting. Clearcutting is practiced only for small areas (lower than 3 ha) in Norway spruce (P. abies) and non-natural forest stands [47,48]. The rest of the forested land is either subjected to conservation management to protect forest health (e.g., phytosanitary felling, trees affected by small local windthrows) or is strictly protected with no management actions allowed. The most common natural disturbances in the area include windthrows and insect attacks [49]. Windthrow events mostly affect conifers in the north sector of the study area [50], particularly Picea forests, with insect outbreaks (Hylobius abietis L., Ips duplicatus Sahlberg and I. typographus L.) being a common secondary disturbance following windthrows [51] and drought [49,52] in conifer forests.

Landis-II Model
Landis-II is a collection of spatially explicit forest landscape models [53] that simulate forest change as a function of succession and disturbances [33]. The landscape in Landis-II is defined as a grid of cells (here 200 × 200 m), each of which belongs to an ecoregion (i.e., areas of homogeneous soil and climate) and can contain multiple species cohorts that can be independently killed by disturbances, competition or age-related mortality, as the succession progresses. Landis-II integrates a number of ecological process models through its modular design. Here, we used the PnET succession extension [54] to underpin tree species establishment, growth, mortality and decomposition. This extension embeds elements of the PnET ecophysiology model of [55] and accounts for competition for available light and water. Biomass growth is the result of a number of processes (e.g., photosynthesis, evapotranspiration) controlled by species ecophysiological parameters (e.g., foliar N concentration, photosynthetic rates) given a number of conditions that include precipitation, temperature and atmospheric CO 2 concentration. Climate change and CO 2 enrichment are interwoven as change in environmental parameters over time, making the PnET extension convenient for climate change modeling [30,56,57].
Landis-II interdependent disturbance extensions were used to model the impacts of harvest [31,58], wind [59] and insect outbreaks [60,61]. The harvest module [62] implements prescriptions in different management areas according to a temporal schedule, management systems and stand characteristics, including resource availability, age structure or stand composition. The biological disturbance agent module [63] implements insect impacts based on pest species preferences and resource availability. The wind module [64], which models windthrow events and operates independently of climate, was used to trigger insect outbreaks.

Landscape Design: Ecoregions and Forest Communities
The initial landscape for Landis-II simulations was built based on the local management plans that were available for the public forests (approximately 60% of the total forest land) [37][38][39][40][41][42][43][44][45][46]. Management plans contained the spatial delimitation of forest stands and information of their species composition and cohort ages, with stands ranging from 0.1 to 50 ha, the legal maximum allowed, with a median of 3.25 ha. Stands were classified into 14 ecoregions of homogeneous climate, soil type and tree species abundances (Supplementary Table S1; Supplementary Figure S1). Ecoregions were ascribed to one of the five forest types according to species composition ( Figure 1). To constrain the number of communities, i.e., cohorts and species combinations that conform to a forest stand, a total of 223 initial communities, comprising from 1 to 4 species and a range of 1 to 17 age cohorts, were identified in the management plans and assigned, based on their frequency, to the corresponding stand and ecoregion in the simulated landscape. The area of private forest lands, for which information was not available, was delimited based on Corine Land Cover 2012 [65]. To ascribe ecoregions and initial communities in private forests, a random forest algorithm [66] Forests 2021, 12, 1434 5 of 18 was trained using information from the management plans and environmental variables (DEM, [67], climate [68], soil properties and classification [34] and distance to infrastructure and populations [69,70]). The overall accuracy of the prediction reached 96.7% over the independent testing dataset.

LANDIS-II Parameterization
The modeled species included Abies alba Mill. (silver fir), Alnus glutinosa (L.) Gaertn. (European alder), Alnus incana (L.) Moench (grey alder), C. betulus (European hornbeam), F. sylvatica (European beech), P. abies (Norway spruce), Q. cerris (Turkey oak), Q. frainetto Ten. (Hungarian oak), Q. petraea (Matt.) Liebl. (sessile oak) and Q. robur L. (pedunculate oak). Species contributing less than 2% to the overall wood volume in the study area were not included in the simulation. Tree species parameters required by the model were measured in the study area [71], compiled from unpublished data or obtained from the local Romanian literature [36,72,73] and other international sources [74,75]. Model parameterization was performed following the recommendations of [29]. The three main background disturbances of the studied area were modeled to increase the accuracy of the simulated landscape and were considered a key element in this study.
The wind disturbance extension required the windthrow return interval and severity, which were modeled/introduced based on national historical records [76]. The average wind speed [68] and the area occupied by conifer forests were used to determine the windthrow local impact on each ecoregion. The mortality of age cohorts was based on Popa [76]. It was assumed that wind event occurrence will not be affected by climate change.
Bark beetle outbreaks were linked to windthrows and were followed by harvesting interventions as per standard practices in the study area. Three main bark beetle species were modeled: I. typographus and I. duplicatus, which target mature Picea stands [49,77], and H. abietis, which targets stands with a high density of A. alba saplings [78]. Pest species preferences (host species and age cohort) and impacts were defined based on expert knowledge, forest health status monitoring [79][80][81] and national Romanian research studies [49,51,52,82].
Harvesting prescriptions were defined following the national Romanian regulatory framework [47,48,83,84] which establishes that logging cycles for the main species in the study area are between 100 and 200 years for beech, Norway spruce, sessile oak and silver fir, 70 to 160 years for other oaks and up to 80 years for broadleaved softwood species, depending on the silvicultural system, forest function, site productivity and other specific situations. The simulated landscape was divided into management areas, composed of multiple stands, that included unmanaged, strictly protected forest (<1%), protected forest with special conservation prescription practices (less than 10 m 3 /ha, approximately 40%), selective logging (20%; with approximately half of the area dedicated to production of mixed beech-conifer, mixed beech-broadleaved and oak forests), shelterwood forests (28%; majority of Fagus, also the remaining mixed beech-conifer, mixed beech-broadleaved and oak forests) and clearcut forests (12%; mostly Norway spruce pure stands). Within each management area, management treatments (stand selection, harvesting periodicity, intensity, targeting species and cycle duration) were implemented based on stand composition and stand age.

Climate Change Scenarios
Three climate change scenarios were defined ( Figure 2) based on the RCP 2.6, 4.5 and 8.5 emission scenarios [85]. For the model spin-up period (1901-2015), we used gridded climate data (precipitation, maximum and minimum temperature) from the Climate Research Unit Time Series (CRU TS) [86] high-resolution dataset. For the simulation period between 2015 and 2100, we used RCP projections from Climatic Data Generator (ClimGen) [87] based on the French National Centre for Meteorological Research CNRM-CN5 climate model [88]. For the simulation period between 2100 and 2140, for which there were no cli-Forests 2021, 12, 1434 6 of 18 mate projections available, climate data were randomly resampled from the last 3 decades, the last climate normal period, of the corresponding CNRM-CN5 series. Climate series were extended to 2140 because of the species longevity and the long logging cycles used in Romanian forestry.
were no climate projections available, climate data were randomly resampled from the last 3 decades, the last climate normal period, of the corresponding CNRM-CN5 series. Climate series were extended to 2140 because of the species longevity and the long logging cycles used in Romanian forestry.
Gridded 0.5° climate data were statistically downscaled following Moreno and Hasenauer [89] using a 30′′ resolution gridded climate dataset from WorldClim 1.4 and 2.0 [68]. Mean annual CO2 concentrations taken from van Vuuren et al. [85] differed among PCR scenarios but were assumed spatially constant. Similar to climate, CO2 concentrations after 2100 were kept constant at the level of the last climate normal to prevent inconsistencies with resampled climate series. Photosynthetically active radiation (PAR), calculated from WorldClim 2.0 solar radiation [68], changed seasonally and spatially but did not differ across climate change scenarios.

Model Output and Validation
We simulated changes in forest biomass and species composition over a period of 125 years (2015-2140) under three climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5), using 22 replicates of each scenario at a spatial resolution of the 200 m cell size. Model outputs included annual live, dead and harvested biomass throughout the simulation period. Forest net productivity was calculated as the change in biomass accounting for the harvesting. Model outputs for the period 2010-2020, when the management plans were released, were validated using Romanian forest yield tables [72,73] and forest management plans of the study area [37][38][39][40][41][42][43][44][45][46]. The simulated initial biomass per hectare ( Figure 3) Gridded 0.5 • climate data were statistically downscaled following Moreno and Hasenauer [89] using a 30 resolution gridded climate dataset from WorldClim 1.4 and 2.0 [68]. Mean annual CO 2 concentrations taken from van Vuuren et al. [85] differed among PCR scenarios but were assumed spatially constant. Similar to climate, CO 2 concentrations after 2100 were kept constant at the level of the last climate normal to prevent inconsistencies with resampled climate series. Photosynthetically active radiation (PAR), calculated from WorldClim 2.0 solar radiation [68], changed seasonally and spatially but did not differ across climate change scenarios.

Model Output and Validation
We simulated changes in forest biomass and species composition over a period of 125 years (2015-2140) under three climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5), using 22 replicates of each scenario at a spatial resolution of the 200 m cell size. Model outputs included annual live, dead and harvested biomass throughout the simulation period. Forest net productivity was calculated as the change in biomass accounting for the harvesting. Model outputs for the period 2010-2020, when the management plans were released, were validated using Romanian forest yield tables [72,73] and forest management plans of the study area [37][38][39][40][41][42][43][44][45][46]. The simulated initial biomass per hectare ( Figure 3) and the extracted timber volume were within the ranges recorded in the management plans for all five forest types. As expected, compared to management, background natural disturbances had minor effects on the current landscape. The simulated windthrow impact was in accordance with volumes estimated from Popa [76] and the area affected by intense windthrows during the period 1986-2016 [79][80][81]. Insect outbreak occurrence and impact were also in accordance with the area affected by intense insect attacks in the study area in the period 1986-2016 [79][80][81], and the harvested volumes were in agreement with the management plans [37][38][39][40][41][42][43][44][45][46]. and the extracted timber volume were within the ranges recorded in the management plans for all five forest types. As expected, compared to management, background natural disturbances had minor effects on the current landscape. The simulated windthrow impact was in accordance with volumes estimated from Popa [76] and the area affected by intense windthrows during the period 1986-2016 [79][80][81]. Insect outbreak occurrence and impact were also in accordance with the area affected by intense insect attacks in the study area in the period 1986-2016 [79][80][81], and the harvested volumes were in agreement with the management plans [37][38][39][40][41][42][43][44][45][46].

Statistical Analysis
Biomass results were spatially aggregated for the five main forest types. Differences across scenarios and forest types were assessed by linear mixed models (LME) [90], with changes over time fit to polynomial functions (up to 5 terms) and simulation replicates as a random factor. Natural disturbance outputs (i.e., mean damaged area, number of killed cohorts and severity) and harvested biomass were analyzed for the entire study area. Shifts in forest species composition at the landscape level were analyzed through the change in species biomass, using constrained redundancy analysis [91] with year and climate change scenario as constraining variables.

Disturbances
Windthrows affected a minor percentage of the forest area (1%), with all of the effect located in the Norway spruce forest type. The affected area significantly increased over time (p ≤ 0.05) such that by 2140, the area damaged was two-fold the area of 2015 ( Figure  4a). A similar trend was observed in the severity and the number of cohorts killed by wind, with values increasing over time, and the number of cohorts killed being significantly lower for the RCP 8.5 scenario. These results are in accordance with the aging of Norway spruce forests, as the susceptibility to windthrows increases with tree age.

Statistical Analysis
Biomass results were spatially aggregated for the five main forest types. Differences across scenarios and forest types were assessed by linear mixed models (LME) [90], with changes over time fit to polynomial functions (up to 5 terms) and simulation replicates as a random factor. Natural disturbance outputs (i.e., mean damaged area, number of killed cohorts and severity) and harvested biomass were analyzed for the entire study area. Shifts in forest species composition at the landscape level were analyzed through the change in species biomass, using constrained redundancy analysis [91] with year and climate change scenario as constraining variables.

Disturbances
Windthrows affected a minor percentage of the forest area (1%), with all of the effect located in the Norway spruce forest type. The affected area significantly increased over time (p ≤ 0.05) such that by 2140, the area damaged was two-fold the area of 2015 (Figure 4a). A similar trend was observed in the severity and the number of cohorts killed by wind, with values increasing over time, and the number of cohorts killed being significantly lower for the RCP 8.5 scenario. These results are in accordance with the aging of Norway spruce forests, as the susceptibility to windthrows increases with tree age.  Insect outbreaks annually affected around 200 ha of the forest area (~1%), mostly in the spruce and mixed Fagus (beech)-conifer forests. The area and the number of cohorts killed by Ips insects, which target mature P. abies trees, increased over time (p < 0.05), also reflecting a change in the forest age structure. This trend was not observed for H. abietis (Figure 4b).
Harvesting occurred mostly in the Fagus and Quercus forest types, where annual means of 300,000 and 100,000 tons were extracted, respectively, with extraction rates of about 10 and 15% every 10 years (Figure 4c). Harvesting in the spruce and mixed beechconifer forest types was mostly restricted to conservation works, with rates of extraction below 2%. Differences between climate change scenarios were noted from 2070, with biomass harvested in RCP 8.5 being greater, followed by RCP 4.5, and RCP 2.6, the lowest among all the scenarios. Given that harvesting prescriptions were constant over time, following the national Romanian regulatory framework, this indicates that differences in forest attributes developed over time among climate change scenarios.

Changes in Biomass across Climate Change Scenarios
Living aboveground biomass ( Figure 5) tended to increase in all forest types and climate change scenarios throughout the simulation period, with fluctuations related to disturbance or harvesting events. Comparatively larger increases in biomass were observed for RPC 8.5, ranging from 21% in the Quercus forest type to 51% in the mixed beech-broadleaved forest type, than for RCP 4.5 (from 2% to 37%). The lowest change in biomass was observed in RCP 2.6, ranging from a decrease in total biomass (−15%) in the Quercus forest type to a 29% increase in the mixed beech-conifer forest type. The effect of climate varied across forest types: mixed beech-broadleaved forests showed the largest differences in the final biomass among scenarios, with RCP 8.5 showing a larger biomass than RCPs 4.5 and 2.6.
The analysis of the biomass at the species level revealed changes in species biomass through the simulation timespan and among climate change scenarios ( Figure 6). The five Insect outbreaks annually affected around 200 ha of the forest area (~1%), mostly in the spruce and mixed Fagus (beech)-conifer forests. The area and the number of cohorts killed by Ips insects, which target mature P. abies trees, increased over time (p < 0.05), also reflecting a change in the forest age structure. This trend was not observed for H. abietis (Figure 4b).
Harvesting occurred mostly in the Fagus and Quercus forest types, where annual means of 300,000 and 100,000 tons were extracted, respectively, with extraction rates of about 10 and 15% every 10 years (Figure 4c). Harvesting in the spruce and mixed beech-conifer forest types was mostly restricted to conservation works, with rates of extraction below 2%. Differences between climate change scenarios were noted from 2070, with biomass harvested in RCP 8.5 being greater, followed by RCP 4.5, and RCP 2.6, the lowest among all the scenarios. Given that harvesting prescriptions were constant over time, following the national Romanian regulatory framework, this indicates that differences in forest attributes developed over time among climate change scenarios.

Changes in Biomass across Climate Change Scenarios
Living aboveground biomass ( Figure 5) tended to increase in all forest types and climate change scenarios throughout the simulation period, with fluctuations related to disturbance or harvesting events. Comparatively larger increases in biomass were observed for RPC 8.5, ranging from 21% in the Quercus forest type to 51% in the mixed beech-broadleaved forest type, than for RCP 4.5 (from 2% to 37%). The lowest change in biomass was observed in RCP 2.6, ranging from a decrease in total biomass (−15%) in the Quercus forest type to a 29% increase in the mixed beech-conifer forest type. The effect of climate varied across forest types: mixed beech-broadleaved forests showed the largest differences in the final biomass among scenarios, with RCP 8.5 showing a larger biomass than RCPs 4.5 and 2.6. main species, P. abies, A. alba, F. sylvatica, Q. petraea and Q. frainetto, showed significant biomass increments during the succession (F. sylvatica up to 2050) (p < 0.05). Of those five, all but F. sylvatica showed significant differences in biomass among the scenarios, with the largest biomass under RCP 8.5 (p < 0.05). Three species, C. betulus, Q. petraea and particularly F. sylvatica, showed short periods with strong reductions in biomass after 2050, which co-occurred in time with low SPEI values, although they tended to recover later.  The analysis of the biomass at the species level revealed changes in species biomass through the simulation timespan and among climate change scenarios ( Figure 6). The five main species, P. abies, A. alba, F. sylvatica, Q. petraea and Q. frainetto, showed significant biomass increments during the succession (F. sylvatica up to 2050) (p < 0.05). Of those five, all but F. sylvatica showed significant differences in biomass among the scenarios, with the largest biomass under RCP 8.5 (p < 0.05). Three species, C. betulus, Q. petraea and particularly F. sylvatica, showed short periods with strong reductions in biomass after 2050, which co-occurred in time with low SPEI values, although they tended to recover later.   The percentage of dead wood biomass relative to the total biomass showed a slight increase in all the forest types over time (Figure 7). Whereas the living biomass of Norway spruce and Quercus spp. differed among scenarios, the RCP scenario did not affect the percentage of dead biomass in these forests, which also showed a trend with very small interdecadal variations. Forest types with participation of F. sylvatica were subjected to strong variations in dead biomass over time, usually within the range 2.5-5.0%, which is also related to fluctuations in living biomass ( Figure 5) and a low 6-month SPEI (Figure 2). The mortality events after drought temporally altered the proportion of live and dead biomass and occurred with a higher intensity in the RCP 8.5 scenario. Figure 6. Mean (dashed line) and 95% confidence interval (shaded area) of the species living biomass (·10 3 kg·ha −1 ) measured over time under the three climate change scenarios: RCP 2.6 (green), RCP 4.5 (orange) and RCP 4.5 (red).
The percentage of dead wood biomass relative to the total biomass showed a slight increase in all the forest types over time (Figure 7). Whereas the living biomass of Norway spruce and Quercus spp. differed among scenarios, the RCP scenario did not affect the percentage of dead biomass in these forests, which also showed a trend with very small interdecadal variations. Forest types with participation of F. sylvatica were subjected to strong variations in dead biomass over time, usually within the range 2.5-5.0%, which is also related to fluctuations in living biomass ( Figure 5) and a low 6-month SPEI ( Figure  2). The mortality events after drought temporally altered the proportion of live and dead biomass and occurred with a higher intensity in the RCP 8.5 scenario.

Productivity
The main forest types showed different aboveground net productivities, with Norway spruce and mixed forest productivity being typically lower than 15 kg/ha, that of beech pure forests being around that value and that of Quercus being usually above it. The aboveground net productivity tended to decrease slightly over time for all forest types. Productivity also differed among climate change scenarios, with the magnitude of the differences depending on forest type (Figure 8). Three forest types, Norway spruce, mixed beech-broadleaved and Quercus, showed significantly (p < 0.05) higher productivity under RCP 8.5 than under the other two scenarios. Beech forest productivity, despite having high temporal variations, showed significant differences between RCP 8.5 and RCP 4.5 (p < 0.02), usually being larger in RCP 8.5. The mixed beech-conifer forest productivity only showed significant differences between the RCP 2.6 and RCP 4.5 scenarios (p = 0.035).

Productivity
The main forest types showed different aboveground net productivities, with Norway spruce and mixed forest productivity being typically lower than 15 kg/ha, that of beech pure forests being around that value and that of Quercus being usually above it. The aboveground net productivity tended to decrease slightly over time for all forest types. Productivity also differed among climate change scenarios, with the magnitude of the differences depending on forest type (Figure 8). Three forest types, Norway spruce, mixed beech-broadleaved and Quercus, showed significantly (p < 0.05) higher productivity under RCP 8.5 than under the other two scenarios. Beech forest productivity, despite having high temporal variations, showed significant differences between RCP 8.5 and RCP 4.5 (p < 0.02), usually being larger in RCP 8.5. The mixed beech-conifer forest productivity only showed significant differences between the RCP 2.6 and RCP 4.5 scenarios (p = 0.035).

Changes in Landscape Species Composition
Forest species composition varied over time and across climate change scenarios (Figure 9). During the first half of the simulations, the landscape species composition was similar among scenarios; in the following decades, the position of samples started to drift with the displacement increasing over time. This resulted in a significant differentiation in the species composition between RCP 8.5 and the other two scenarios. As shown by the ordination, and in accordance with the biomass analysis, RCP 8.5 was characterized by a greater biomass of Q. frainetto, P. abies and A. alba.

Changes in Landscape Species Composition
Forest species composition varied over time and across climate change scenarios (Figure 9). During the first half of the simulations, the landscape species composition was similar among scenarios; in the following decades, the position of samples started to drift with the displacement increasing over time. This resulted in a significant differentiation in the species composition between RCP 8.5 and the other two scenarios. As shown by the ordination, and in accordance with the biomass analysis, RCP 8.5 was characterized by a greater biomass of Q. frainetto, P. abies and A. alba.   Fagus-conifers, Fagus, Fagus-broadleaved, Quercus) for the forecasts for every decade and climate change scenario: RCP 2.6 (green), RCP 4.5 (orange) and RCP 4.5 (red).

Changes in Landscape Species Composition
Forest species composition varied over time and across climate change scenarios (Figure 9). During the first half of the simulations, the landscape species composition was similar among scenarios; in the following decades, the position of samples started to drift with the displacement increasing over time. This resulted in a significant differentiation in the species composition between RCP 8.5 and the other two scenarios. As shown by the ordination, and in accordance with the biomass analysis, RCP 8.5 was characterized by a greater biomass of Q. frainetto, P. abies and A. alba.

Discussion
The simulation model implemented for temperate forests in the Romanian Southern Carpathians showed that climate change produces changes in forest biomass, productivity and species abundance. Such changes were particularly noticeable under RCP 8.5, the most extreme climate change scenario, and rendered an overall increase in the carbon carrying capacity of the studied Carpathian forests.

Forest Biomass and Productivity
The results of our Landis-II model adapted to Southern Carpathian forests showed an overall increase in forest biomass over time, with significantly greater biomass accumulation in the RCP 8.5 scenario. This is in line with Hubau et al. [10], who found that contrary to some tropical forests, temperate forests maintain a certain capacity to keep stocking carbon. The increased forest biomass accumulation and productivity predicted for the study area with intense climate change point to the raised CO 2 concentration [2,5] and warmer winters [3] as being the main drivers in the Carpathians, where precipitations do not change significantly for most of the area [23]. However, our simulations also indicate that the initial increase in biomass is limited, likely due to increasing drought and the stabilization of climate and atmospheric CO 2 concentrations. This trend suggests that the equilibrium already observed in some tropical forests [10] might occur in temperate forests somewhere at the end of the current century. Indeed, other authors have reached similar conclusions and timing, pointing directly to the proliferation of bark beetles as one of the main natural disturbances limiting carbon sequestration associated with climate change [92]. In the Southern Carpathians, climate change conditions overcompensate the negative effects of the disturbances at the beginning of the period but may not be able to do that after 2040, for most of the forests.
The main forest types studied responded differently to climate change, according to the ecophysiology of the species. Norway spruce and Quercus spp. forests had low productivity but higher stability, while mixed beech-conifer forests, beech forests and mixed beech-broadleaved forests were more prone to biomass changes, particularly after 2050. Indeed, low-and medium-altitude forests in the Southern Carpathians are very likely to be more affected by drought, while the precipitation in mountain tops, often covered by Norway spruce forests, might increase slightly [23]. Bouriaud el al. [16] also predicted forest biomass increments for the Romanian Eastern Carpathians but found a decreasing trend for conifer forests. This discrepancy among Bouriaud et al.'s [16] results and ours is likely due to three causes: (a) climate change effects on forests vary in type and intensity at the global, regional and local levels [15,23,93], as do the local conditions and climate change projections of the study areas in the Eastern and Southern Carpathians; (b) the modeling approach, with PnET and Landis-II models also accounting for complex ecophysiological responses to photosynthetically active radiation, atmospheric CO 2 concentration, etc.; and (c) the climate models used [16,[94][95][96][97]. Indeed, the CNRM-CN5 model [87,88] used here is one of the least limiting climate models regarding water availability in the study area.
The predicted increase in the impact of the main natural disturbances in the Southern Carpathians: windthrows and bark beetle attacks, particularly Ips sp., is in agreement with other studies [6,98]. Even though some variations may be caused by harvest and aging forests, the periods with low productivity are likely associated with drought, since the SPEI value tends to decrease within the same periods. Drought events also occur simultaneously with raised mortality and deadwood biomass increments. This is the case of Norway spruce forests and beech forests after 2040. Natural disturbances induce organic carbon release, and hence it is common that disturbance regimes and ecosystem resilience determine forests' carbon sequestration capacity [10,[99][100][101][102]. A net increase in temperate forest biomass caused by climate change has been forecasted for areas where disturbances and extreme climate events have a limited impact [54,103]. These disturbances usually damage conifer forests [6], and aging forests will likely contribute to it [52]. In our model, H. abietis had a minor impact, as a consequence of the low biomass of the cohorts they target, and because the selection harvesting management implemented in mixed beech and silver fir (A. alba) stands prevents H. abietis outbreaks [104]. The change in wood harvest observed here was also predicted by Bouriaud et al. [16], Ciceu et al. [105] and Chivulescu et al. [33] and concentrated in the intermediate-and low-altitude forests, having consequences for the forest structure. This increment in harvested timber is directly associated with both forest aging and the application of logging cycles [78]. As a result of all previous processes, around the year 2050, forests are expected to reach the carbon sequestration limit. Such an idea is supported by Bouriaud et al.'s [16] and Hubau et al.'s [10] studies. Nabuurs et al. [19] also found evidence of carbon sink saturation in European forests.

Forest Composition
The Landis-II projections predicted shifts in forest species composition and abundance under the RCP scenarios in Romanian Southern Carpathian temperate forests, with the five main species showing significant biomass increments over time (F. sylvatica only up to 2040). Indeed, the reduction in the cohorts killed by windthrows under RCP 8.5 compared to RCP 4.5 and RCP 2.6 is likely due to the success of A. alba and the intense modification of the forest structure under RCP 8.5. However, only F. sylvatica showed differences among the RCP scenarios, with the highest biomass accumulation under RCP 8.5. Some species, particularly F. sylvatica, showed short periods with strong reductions in the biomass after 2050 that co-occurred, and thus are likely associated, with drought conditions. Climate change may alter local environmental conditions, making them no longer suitable for a given species and altering the interspecific competition [18,96], ultimately resulting in a change in species abundance [1]. For instance, Musselman and Fox [14] indicated that drought-tolerant species are expected to succeed at the global scale under the new conditions, and this is likely going to occur for A. alba in the Southern Carpathians, as predicted here.
Previous studies have indicated that the intensity of the change in climate is proportional to the magnitude of the impact on forests [16,33]. The intensity of such a change in the Southern Carpathians is expected to be relatively small, as the climate conditions do not strongly depart from the current conditions under the analyzed RCPs. However, some forest types are more vulnerable to changes than others [18,96]. Norway spruce forests at the upper limit will gain from the temperature rise [15], and, as suggested by Landis-II models, at the lower limit, they will also benefit from F. sylvatica's vulnerability to drought. At lower altitudinal belts, conifers coexist with beech. Beech and conifer mixed forests are more successful than pure beech stands [25,26], but Norway spruce grows more and remains relatively unaffected compared to F. sylvatica during drought [106][107][108]. A. alba is also less susceptible to warmer and drier conditions than F. sylvatica [25]. Thus, at the stand level, the association between F. sylvatica and A. alba is not only advantageous for enduring drought [25,26,106,107] but also for sunlight use [104,109]. Our results for the Southern Carpathians suggest that climate change will contribute to A. alba encroachment in F. sylvatica pure forest stands. Indeed, there is strong evidence that beech forests have already been affected by climate change-induced drought in recent decades [20]. Moreover, Broadmeadow et al. [4] recently highlighted beech's vulnerability in temperate forests across Europe.
Conifer and beech forests aside, Quercus spp.-dominated forests will partially benefit from the temperature rise and F. sylvatica decline. Quercus species are relatively drought resistant. According to the Landis-II simulations, and as also supported by previous studies [110], the most thermophile and drought-resistant species in the Southern Carpathians such as Q. frainetto [36] will benefit from climate change and increase their biomass during succession, while Q. petraea will eventually suffer the impact of drought events.

Conclusions
Here, we simulated the changes in the composition and structure of the temperate forests of the Southern Carpathians under three climate change scenarios. Our results indicate that climate change may contribute, overall, to increasing temperate forest productivity and biomass, mainly due to the combination of increasing temperatures and thus extended growing periods in the uplands and only a mild reduction in rainfall across the landscape. Fluctuations in this trend associated with strong biomass reductions and temporary deadwood rise were related to drought periods. However, it is likely that increased drought events may eventually counteract such positive trends.
Climate change selectively affected areas and species, thus contributing to changes in species abundance. P. abies benefitted from warmer conditions, whereas forests dominated by F. sylvatica were vulnerable to climate change, with drought periods associated with large mortality events. The A. alba abundance increased under the new climate conditions at the expense of Fagus's decline. Our results suggest that under the upcoming climate, mixed F. sylvatica-conifer (namely Norway spruce and silver fir) stands will have greater resistance and resilience than those of pure F. sylvatica stands, a finding that will help in guiding future forest management practices.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f12111434/s1, Figure S1: Distribution of the ecoregions within the study area. Ecoregions represent areas of homogeneous climate, soil and forest types, Table S1: Characteristics of the defined ecoregions within each forest type: main species composition, elevation range, soil attributes and main natural disturbances affecting the ecoregion. Disturbances: It: I. typografus; Id: I. duplicatus; Ha: H. abietis [34,35,111,112].

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.