Can mixed forests sequester more CO2 than pure forests in future climate scenarios? A case study of Pinus sylvestris combinations in Spain

Adapting forests to climate change is a critical issue for forest management. It requires an understanding of climate effects on forest systems and the ability to forecast how these effects may change over time. We used Spanish Second National Forest Inventory data and the SIMANFOR platform to simulate the evolution of CO2 stock (CO2 Mg · ha−1) and accumulation rates (CO2 Mg · ha−1 · year−1) for the 2000–2100 period in pure and mixed stands managed under different Shared Socioeconomic Pathways (SSPs) in Spain. We hypothesized that (1) the more optimistic climate scenarios (SSP1 >  > SSP5) would have higher CO2 stock and accumulation rates; (2) mixed stands would have higher CO2 stock and accumulation rates than pure stands; and (3) the behavior of both variables would vary based on forest composition (conifer–conifer vs. conifer–broadleaf). We focused on Pinus sylvestris L., and its main mixtures with Pinus nigra, Pinus pinaster, Fagus sylvatica and Quercus pyrenaica. The SSP scenarios had correlating CO2 stock values in which SSP1 > SSP2 > SSP3 > SSP5, ranging from the most optimistic (SSP1) to the most pessimistic (SSP5). Though pure stands had higher CO2 stock at the beginning, differences with regard to mixed stands were drastically reduced at the end of the simulation period. We also found an increase in the aboveground CO2 proportion compared to belowground in conifer–broadleaf mixtures, while the opposite trend occurred in conifer–conifer mixtures. Overall CO2 accumulation rates decreased significantly from the beginning to the end of the simulation period, but our results indicated that this decline would be less drastic in mixed stands than in pure ones. At the end of the simulation period, CO2 accumulation rates were higher in mixed stands than in pure stands for all mixtures, fractions (aboveground and belowground) and SSPs. Knowing the evolution of mixed forests in different climate scenarios is relevant for developing useful silvicultural guidelines in the Mediterranean region and optimizing forestry adaptation strategies. Better understanding can also inform the design of management measures for transitioning from pure stands to more resource efficient, resistant and resilient mixed stands, in efforts to reduce forest vulnerability in the face of climate change. This work highlights the importance and benefits of mixed stands in terms of CO2 accumulation, stand productivity and species diversity.


Introduction
Global warming is altering environmental conditions around the word. In many landscapes, droughts are expected to increase in frequency and intensity, and new extreme climate events are predicted (Allen 2019). Changes in forest productivity and species distribution have also been reported (Fernandez-de-Una et al. 2015). In light of such prospects, identifying abiotic factors that can limit growth and quantifying local-scale productivity changes make it possible to develop site productivity gradients (Toïgo et al. 2015) as a 1 3 starting point for creating management guidelines. Previous studies have indicated that climate change effects are expected to be harsher in the Mediterranean region (Astigarraga et al. 2020;Martín-Benito et al. 2010). However, Mediterranean species have a wide range of water-use strategies that could play a key role in coping with more extreme conditions. Similarly, complementary mechanisms for water use in mixed forests can enhance forest resilience to extreme drought episodes (Muñoz-Gálvez et al. 2021). In recent decades, these adaptive capabilities have attracted the attention of researchers and managers (Condés et al. 2018;Del Rio et al. 2016), some of whom have proposed mixed forests as an option for mitigating the effects of climate change (Toïgo et al. 2015). Research on mixed forests vs monocultures has found that mixed forests provide greater stability (Muñoz-Gálvez et al. 2021) and even higher productivity (Forrester 2014;Pretzsch 2009; Pretzsch and Schütze 2016) due complementary (Toïgo et al. 2015;Pretzsch and Schütze 2016;Del Río et al. 2017;Riofrío et al. 2017a, b) and more efficient use of resources (Pretzsch 2014;Pretzsch and Schütze 2016;Riofrío et al. 2017a, b). In addition, mixed stands have been reported to have greater resilience and resistance to biotic and abiotic disturbances compared to pure stands (Del Río et al. 2009Río et al. , 2017Pardos et al. 2021).
Scots pine (Pinus sylvestris) is the most widespread Pinus species, covering 10,000 km longitudinally from Spain to Russia (Krakau et al. 2013). There are approximately 1 000 000 ha of pure stands in Spain alone (del Río et al. 2009). This species has great commercial and ecological importance and has been studied extensively in Europe (Krakau et al. 2013). It grows throughout the continent Pretzsch et al. 2020) and into the Mediterranean area, which marks the limit of its natural distribution (Fernandez-de-Una et al. 2015). Although Scots pine often grows in natural and forested pure stands , it can be found in mixed stands with diverse coniferous and broadleaf species. In this study, we focused on four different mixtures (Pinus sylvestris -Pinus nigra, Pinus sylvestris -Pinus pinaster, Pinus sylvestris -Fagus sylvatica and Pinus sylvestris -Quercus pyrenaica) to simulate CO 2 stock and accumulation rates under different climate scenarios.
Pinus sylvestris -Pinus nigra stands usually appear closer to the southern distribution of Scots pine, which coincides with the ecological optimum for black pine (Barbéro et al. 1998). Various growth and yield models have been fitted for Scots pine -Black pine mixed stands in Spain (Trasobares et al. 2004a, b, c), some of which even included non-wood products such as fungi (De Aragón et al. 2007;Herrero et al. 2019;Palahí et al. 2009). However, little is known about species complementarity and interactions.
Pinus sylvestris -Pinus pinaster mixing effects have also been studied. Crown complementarity and vertical stratification in the canopy make this mixture more productive due to higher light interception and use compared to pure stands of either species (Riofrío et al. 2017a, b). This mix was also found to have a positive effect on understory richness and tree regeneration compared to monospecific stands (López-Marcos et al. 2020).
Pinus sylvestris -Fagus sylvatica stands are probably the most extensively studied mixed forests in Europe. The vast distribution of these mixed stands is due to species niche complementarity (Del Río et al. 2017;Pretzsch and Schütze 2016) in terms of space occupancy efficiency (Pretzsch andSchütze 2009), light tolerance, light use (Del Río et al. 2016;Pretzsch et al. 2015a, b), structural and vertical heterogeneity (Del Río et al. 2017;Pretzsch and Schütze 2009) and root systems Pretzsch et al. 2015a, b;Yeste et al. 2021). This generally translates into higher productivity and greater water use efficiency compared to pure stands, of either species though stand structure and climate conditions influence outcomes in complex ways (Condés et al. 2013(Condés et al. , 2018Del Río et al. 2017;González-de-Andrés et al., 2017Pretzsch et al. 2015a, b).
Pinus sylvestris -Quercus pyrenaica stands are interesting because of their complementarity and long-term stability in average climate conditions. Pyrenean oak shows low resistance to drought with high recovery rates, while Scots pine has the opposite behavior (Muñoz-Gálvez et al. 2021). Differences in shade tolerance, leaf habits and root depth (Martín-Gómez et al. 2017) lead to higher resource use efficiency in these mixed stands compared to monocultures of either species (Forrester 2014).
Given the importance of Scots pine across Europe and the fact that the Mediterranean region is considered the ecological limit of its distribution (Fernández-de-Una et al. 2015), mixed stands of Scots pine with the aforementioned species hold great interest for the development of Spanish management guidelines. Most studies on the benefits of these four species in mixed forests with Pinus sylvestris have focused on the historical evolution of existing stands, using past experience to guide future forest management (Aguirre et al. 2019;Pretzsch et al. 2015a, b;Steckel et al. 2019). While this is quite reasonable, it misses the point that future climate conditions are likely to be different to those of the last 100 years. In Spain, global warming implies that the frequency and intensity of drought events are expected to increase (Kjellström 2004). Thus, empirical guidance based on past data may prove inadequate for future realities.
Forest management techniques can be adapted to respond to the need for mixed forests as a promising way to reduce forest drought stress (Pretzsch et al. 2015a, b;Steckel et al. 2019).
This study was designed to simulate the evolution of CO 2 stock and accumulation rates in pure and mixed stands under different climate change scenarios for the period spanning from the years 2000 to 2100.
We hypothesized that: (1) The more optimistic climate scenarios (SSP1 > > SSP5) would have higher CO 2 stock and accumulation rates; (2) mixed stands would have higher CO 2 stock and accumulation rates than pure stands; and (3) the behavior of both variables would vary according to forest composition (conifer-conifer vs conifer-broadleaf).

Material and methods
The simulations in this study were developed using data from the Spanish Second National Forest Inventory (SNFI2), which is the most extensive spatio-temporal forestry database in Spain, along with new climate-sensitive models and the SIMANFOR platform. In addition to these innovative and robust resources, WordClim was used, a well-recognized climate database that covers our study area.

Data
Forestry data were extracted from Spanish Second National Forest Inventory (SNFI2) plots. These consist of four concentric circles with radii of 5,10,15 and 25 m, in each of which multiple tree-level variables were recorded for all trees over 7.5, 12.5, 22.5 and 42.5 cm diameter at breast height (1.3 m), respectively Herrero and Bravo 2012). Expansion factors were used to estimate stand variables from individual tree variables such as density (N), quadratic mean diameter (Dq), dominant height (Ho), and aboveground (Wa) and belowground (Wr) biomass. Although this study focuses on mixed forest stands, we selected both pure and mixed SNFI2 plots to look at how mixing effects influence tree and stand dynamics. SNFI2 plots (composed of two main species) were considered mixed when the combined proportion of both species (based on basal area) was greater than 90% and the proportion of either species was over 15%. Pure SNFI2 plots (basal area of the main species 90% or higher) found close to the mixed plots were also selected to create triplets for comparison, consisting of the mixed plot and one pure plot of each of its constituent species. Table 1 presents the main tree and stand variables for the selected SNFI2 pure and mixed plots, and Fig. 1 shows their location. The selection criteria prioritized finding mixed plots for each of the four mixtures in the study that also had some pure stands of the constituent species nearby for comparison. Accordingly, useful plots were located in areas where climate and soil conditions allowed for the presence of both species.

CO 2 stock and accumulation rate simulation in Spanish mixed forests
In this study, using the SIMANFOR platform ), CO 2 stock was estimated for five-year periods from 2000 to 2100 under four different SSPs (SSP1, SSP2, SSP3 and SSP5). SIMANFOR was developed to simulate sustainable forest management alternatives. It allows users to compare different silvicultural scenarios and choose one that best fits their management objectives. From forest inventory and silvicultural scenario inputs, inventory data are processed using equations for the selected model. Successive steps provide the user with output information about individual tree and stand evolution at every stage of the planned scenario. Until now, the SIMANFOR simulator has only worked for pure stand models without climate influence. For this study, new tree growth and allometric functions were added to the SIMANFOR platform, making it possible to simulate dynamics in mixed forest stands under various climate conditions. Tree growth was characterized by tree basal area increment (BAI) and estimated using the distance-independent climate-dependent mixed models fitted by Rodríguez-de-Prado et al. (2022a), which are presented in Table 2.
Diameter at breast height was then derived from BAI estimates and used to predict total tree height using h-d models (Rodríguez-de-Prado et al. 2022b), as shown in Table 3.
We controlled for mortality based on the maximum stand carrying capacity (SDI max ), which defines the maximum number of trees per hectare where natural mortality occurs in a forest stand. At each iteration, the Stand Density Index (SDI) (Reineke 1933) and SDI max were estimated using the climate-dependent MSDR models found in Rodríguez- de-Prado et al. (2020) for each species in each SNFI2 sample plot. Thus, while SDI max ≥ SDI, mortality does not occur in a specific plot for a given species, when SDI max < SDI, the expansion factor for all the trees of that species decreases by 2%. As a consequence, the stand density value (trees · ha −1 ) of the species decreases by 2% and the simulation assumes that change going forward. Biomass compartments (root, stem, leaves, crown) were calculated for each tree using equations from Ruiz-Peinado et al. (2011. Biomass estimates were grouped into above-and belowground biomass and then transformed into C values using the corresponding conversion factors for each species: Fagus sylvatica (0.486), Pinus nigra (0.464), Pinus pinaster (0.468), Pinus sylvestris (0.459) and Quercus pyrenaica (0.457). Finally, the C values were multiplied by 3.67 to obtain CO 2 values. The simulations assumed null incorporated forest mass or regeneration processes. All the SIMANFOR simulations were done using the supercomputing services of a regional Scientific Computing Center known as SCAYLE (SCAYLE 2021). Its higher computational capacity compared to desktop and web versions increased the velocity of the simulation process. After that, all the generated outputs (at tree and plot level) were combined, re-structured and analyzed using R language programming (R Core Team 2020).

Shared socioeconomic pathways (SSPs)
Four different Shared Socioeconomic Pathways scenarios (SSPs) were considered in this paper. Shared Socioeconomic Pathways are scenarios describing how the world would look in the absence of climate policy. They allow researchers to examine barriers and opportunities related to climate mitigation and adaptation by combining possible future scenarios with mitigation targets (Meinshausen et al. 2020;Riahi et al. 2017). SSPs are based on five narratives (Fig. 2) representing distinct future socioeconomic projections and political environments that span the range of plausible futures: SSP1 and SSP5 anticipate reasonably positive patterns in human development, including "significant investments in education and health, robust economic growth, and well-functioning institutions." However, they differ in that SSP5 assumes this would be powered by an energyintensive, fossil-fuel economy, while SSP1 assumes a gradual transition toward sustainable practices. SSP3 and SSP4 are more negative about the economic and social growth potential of their nations; they reflect low investment in education and health in poorer countries along Table 3 Tree height-diameter equations for individual and paired species h total tree height (m), d diameter at breast height (cm), do dominant diameter (cm), Ho dominant height (m), dq quadratic mean diameter (cm), dq i quadratic mean diameter of the species i (cm), BAL basal area of larger trees (m 2 ·ha −1 ), Pi species proportion by area, M De Martonne Aridity Index (mm·°C −1 ), β 0 and β 1 model parameters to be estimated, RMSE root mean square error

Species composition Species Equation RMSE
Pinus sylvestris -Pinus nigra Pinus sylvestris

Pinus sylvestris -Fagus sylvatica
Pinus sylvestris with rapidly rising population and growth disparities. SSP2 is a "middle of the road" scenario in which historical growth trends are maintained in the twenty-first century. Data from SSP4 were not available in WorldClim, so this scenario was not analyzed in our study.

Pinus sylvestris -Quercus pyrenaica Pinus sylvestris
The De Martonne aridity index (De Martonne 1926) calculated as P/(T + 10) (where P is the total annual precipitation in mm, and T the mean annual temperature in °C) is a climate index commonly still used to describe the aridity gradient or drought in a given area (Aguirre et al. 2018;Bielak et al. 2014;Condés et al. 2017). In this study, the De Martonne Index was one of the variables selected to predict mortality in the SDI and SDI max models (Rodríguez-de-Prado et al. 2021) as well as tree growth for BAI productivity and h-d models (Tables 2 and  3). The climate scenarios described above were used to predict four different De Martonne aridity indexes in each SFNI2 plot. Each plot was simulated four times (one per SSP scenario), and four De Martonne aridity indexes were calculated for each plot, based on the climate predictions for four time periods: 2000-2020, 2020-2040, 2040-2060 and 2060-2100.

Forest CO 2 stock simulation for Pinus sylvestris mixtures during the 2000-2100 period
Simulated CO 2 stock with different Pinus sylvestris mixtures and SSPs for the 2000-2100 period is shown in Supplementary Table 4. We found a common SSP1 > SSP2 > SSP3 > SSP5 trend with higher stock values in the most optimistic scenario (SSP1) and lower values in the most pessimistic scenario (SSP5). The results also show that CO 2 stock patterns evolved differently in pure and mixed stands. Early in the simulation period, mixed stands generally had less stock capacity than their respective pure stands. However, our findings indicate that the differences reduced drastically as time progressed and that by the end of the simulation period mixed stands had greater CO 2 stock capacity than pure stands.
Based on SSP2 as an intermediate future scenario, Pinus sylvestris-Fagus sylvatica mixed stands had less stock capacity than Pinus sylvestris (1.0%) and Fagus sylvatica (57.5%) pure stands in the year 2000 (Fig. 3). By the year 2100, these differences had changed to 3.6% and 9.8%, respectively. Similarly, Pinus sylvestris and Pinus nigra pure stands both had ~ 17.5% higher stock rates than their mixture early on. By the end of the simulation period, stock in the mixed stands was 6.4% higher than in Pinus sylvestris pure stands and the difference with respect to Pinus nigra monocultures had decreased by ~ 10.0%. In the Pinus sylvestris-Pinus pinaster mixture, the mixed stands initially presented greater stock capacity (11.0%) than Pinus pinaster pure stands but less than Pinus sylvestris pure stands (27.8%). By the year 2100, the mixed stands had greater stock capacity (15.9% and 8.1%) than Pinus pinaster and Pinus sylvestris pure stands, respectively. A similar trend was found in Pinus sylvestris-Quercus pyrenaica mixtures, with initially higher CO 2 stock in the mixed stands compared to Quercus pyrenaica pure stands (9.7%) and lower CO 2 stock compared to Pinus sylvestris pure stands (47.5%). Differences with the latter reduced in successive timespans to the point that CO 2 stock in mixed stands surpassed that of pure stands. At the end of the simulation period, the mixed stands had 25.7% and 11.7% more stock capacity than Quercus pyrenaica and Pinus sylvestris pure stands, respectively.
Aboveground and belowground CO 2 were also simulated in our study (Supplementary Table 4), and several trends were observed based on species traits in the mixtures analyzed. In conifer-broadleaf mixtures, we found that the aboveground proportion of stored CO 2 tended to increase in successive simulation periods in both pure and mixed stands. For these mixtures, aboveground biomass increased by 1.1% in Pinus sylvestris-Fagus sylvatica mixed stands and by 2.2% in Pinus sylvestris-Quercus pyrenaica mixed stands during the study period.
The opposite trend occurred for conifer-conifer mixtures, both of which experienced a small but steady decrease in aboveground biomass proportion between 2000 and 2100 in all SSPs. However, differences were found between pure and mixed stands of Pinus sylvestris-Pinus nigra and Pinus sylvestris-Pinus pinaster. In the first, the aboveground stored CO 2 proportion for Pinus sylvestris pure stands decreased over time (2.0%), while Pinus nigra pure stands increased by the same amount. In the second, the aboveground stored CO 2 proportion for Pinus pinaster pure stands decreased notably (3.7%) but increased for Pinus sylvestris pure stands as time progressed.
The results presented in this section are also available in terms of biomass (Mg · ha −1 ) in Supplementary Table 5.

Forest CO 2 accumulation rate simulation for the 2000-2100 period
Based on simulated CO 2 stock values, accumulation rate in terms of megagrams of CO 2 and biomass (above-and belowground) per hectare and year was determined for successive periods from 2000 to 2100 (Supplementary Tables 6 and 7).
In the initial period, total biomass accumulation rates for Pinus sylvestris in mixed stands ranged from 4.0-6.6 Mg · ha −1 · yr −1 . Differences were observed between pure stands and mixtures. Pinus sylvestris -Pinus nigra had similar accumulation rates (4-4.5 Mg · ha −1 · yr −1 ) for mixed and pure stands. Mixed stands of Pinus sylvestris-Fagus sylvatica presented higher accumulation rates than Pinus sylvestris pure stands (23.9%), but slightly lower rates (5.4%) than Fagus sylvatica pure stands. However, Pinus sylvestris mixtures with Pinus pinaster or Quercus pyrenaica presented higher accumulation rates than their respective pure stands (~ 2 Mg · ha −1 · yr −1 ). Similar patterns were found for the aboveground and belowground fractions also, with two exceptions: The accumulation rate in Pinus sylvestris pure stands was initially lower than in Fagus sylvatica and their mixed stands in the aboveground fraction (64.1% and 46.2%, respectively), but the belowground accumulation rate was higher (25.0% and 3.8%, respectively) and Pinus sylvestris pure stands showed slightly higher aboveground accumulation rates than Pinus nigra pure stands and their mixed stands (3.4% and 5.3%, respectively), but the opposite occurred in belowground accumulation rates (68.7% and 20.0%, respectively).
For all analyzed mixtures, accumulation rates for both pure and mixed stands decreased significantly from the beginning to the end of the simulation period under all SSPs included in the study (Supplementary Table 6). Our results indicate that these reductions would be less drastic in mixed stands than in pure ones. Here, it is important to highlight that accumulation rates in all mixed stands, fractions (aboveground and belowground) and SSPs were higher than in pure stands at the end of the simulation period.
Using SSP2 as the "middle-of-the-road" future scenario, differences in accumulation rates between 2000 and 2100 ranged from 30 to 50% in mixed stands and from 50 to 80% in pure stands. Conifer-conifer mixed stands experienced a similar accumulation rate decrease of ~ 50%. In conifer-broadleaf mixtures, accumulation rates decreased by 36.3% in Pinus sylvestris-Fagus sylvatica and by nearly 60.0% in Pinus sylvestris-Quercus pyrenaica mixtures. Our results for CO 2 stock simulations (Supplementary Table 4) were also consistent with the SSP1 > SSP2 > SSP3 > SSP5 pattern of higher accumulation rates under the most optimistic scenario (SSP1) and lower values under the most pessimistic scenario (SSP5).
In comparing accumulation rate decreases between SSP5 and SSP1, we found that in the Pinus sylvestris-Fagus sylvatica mixture, the Pinus sylvestris accumulation rate decreased by 6.0% and Fagus sylvatica by 7.9%, compared to a 19.1% reduction in their mixed stand. In the Pinus sylvestris-Quercus pyrenaica group, the accumulation rate decreased by 3.5% in Pinus sylvestris and by 2.3% in Fagus sylvatica, compared to a 3.4% reduction in their mixed stand. Figure 4 shows the simulated accumulation rates of the various mixtures and forest types (pure and mixed) for the 2000-2100 period. Though conifer-broadleaf mixed stand species had different initial growth trends, the outcomes were quite similar. Fagus sylvatica pure stand accumulation rates for CO 2 were initially higher than Pinus sylvestris and Pinus-Fagus mixed stands for aboveground production (35.3% and 9.1%, respectively) and lower for belowground production (40.9% and 31.2%, respectively). Nevertheless, their accumulation rates decreased faster than the other cases and fell below Pinus sylvestris and mixed stand values for above-and belowground CO 2 production at around 2050 and 2090, respectively. Additionally, mixed stand production was higher than Pinus sylvestris pure stands during the entire study period and accumulation rate reduction over time had a shallower slope. At the end of the simulation (using SSP2 as reference), total CO 2 accumulation rate reductions in mixed stands were lower (36.1%) than in pure Fagus sylvatica (80.0%) or Pinus sylvestris (54.8%) stands. The initial accumulation rate in Pinus sylvestris pure stands was higher (22.1%) than in Quercus pyrenaica monocultures, but by the end of the study period, Quercus pyrenaica pure stands were more productive (41.4%) than Pinus sylvestris stands, as reflected in both the above-(38.1%) and belowground (50.0%) accumulation rates. Nonetheless, Pinus-Quercus mixed-stand accumulation rates were higher than those of their pure stands throughout the studied period, with dramatic differences of 56.7% and 30.0% for aboveground and 61.0% and 21.9% for belowground accumulation rates, respectively, by the year 2100.
The CO 2 accumulation rates in conifer-conifer mixtures led to similar stock results. Pinus nigra had lower accumulation rate reduction (60.8%) than Pinus sylvestris (72.3%) in pure stands, except for the aboveground accumulation rate during the first period. The accumulation rate for Pinus nigra pure stands was 8.5% higher than in mixed stands during the early periods, but mixed stands eventually surpassed them by as much as 15.7% (in the Fagus sylvatica mixture, aboveground fraction). In the other conifer mixture, Pinus pinaster pure stands were more productive than Pinus sylvestris pure stands for aboveground, belowground and total production from the beginning (21.1% of total production) to the end (22.2% of total production) of the study period. However, mixed-stand accumulation rates were always higher than those of the respective pure stands, and the accumulation rate reduction over time was also smaller (50.2% CO 2 accumulation rate reduction). Comparatively, although the initial accumulation rates were different in the two mixtures, the trends described in both cases were quite similar and the aboveground accumulation rate decreased remarkably more than the belowground rate.
While in some cases, initial production was higher in pure stands than in their respective mixtures, the slope for all mixed stands tended to be shallower compared to pure stands (49.2% and 50.1% growth reduction, respectively). Mixedstand above-and belowground productivity surpassed that of the corresponding pure stands at the end of the study period. Additionally, broadleaf-conifer mixtures showed more drastic reductions in belowground (35.6% with Fagus sylvatica and 56.5% with Quercus pyrenaica) than aboveground (39.6% with Fagus sylvatica and 64.4% with Quercus pyrenaica) CO 2 accumulation rate. For conifer-conifer mixtures, however, the opposite occurred in both belowground (49.4% for Pinus nigra and 51.1% for Pinus pinaster) and aboveground (46.5% for Pinus nigra and 48.8% for Pinus pinaster) results. Differences in total CO 2 accumulation rates between the extreme SSPs for mixtures ranged from 1.54% to 10.73%. The latter number corresponds to the mixture with Fagus sylvatica, which proved to be the most extreme case in this study.

Discussion
This study analyzed the simulated evolution of CO 2 stock and accumulation rates in Mediterranean pure and mixed tree stands in four climate scenarios for time period of the year 2000 to 2100. Results showed consistent differences among climate scenarios in CO 2 stock and accumulation rates in the forests studied. Though species mixtures were found to influence stand behavior, the results revealed a common trend of higher CO 2 stock and accumulation rates in mixed vs pure stands.
Our findings also indicated higher CO 2 stock and accumulation rates in the most optimistic climate scenario (SSP1) compared to the most pessimistic one (SSP5). Contrary findings were reported previously for Swedish boreal forests by Poudel et al. (2011), with higher production in the climate change scenarios than under control conditions. Given the initial local climate conditions in those forests, results from that study may support the idea that climate change effects could be favorable there due to low or null rainfall reduction and increased temperatures. In Mediterranean forests, however, the results would indicate a situation of increasing drought stress that decreases CO 2 accumulation rates and subsequent stock capacity. Looking at other factors, Morán-Ordóñez et al. (2020) simulated future carbon stock capacity and other ecosystem services in a Mediterranean forest. They concluded that management policies would be more decisive than climate in providing ecosystem services. In the present work, not applying silviculture, the CO2 accumulation rate decreased on each period, consistently in all SSPs, due to the increase on climate conditions severity, as it was found by (Ma et al. 2014;Steenberg et al. 2011).
Our findings showed that CO 2 stock was higher in pure stands at the start of the study period, but the differences diminished over time and mixed stands even surpassed them in the end. Although species mixtures were found to influence stand responses, a clear tendency toward higher CO 2 stocks in mixed stands emerged over long simulation periods, which is consistent with previous studies, Dai et al. (2016) reported similar findings for subtropical coniferous plantations in China. Muñoz-Gálvez et al. (2021) and Pardos et al. (2021), associated higher CO 2 storage in mixed stands with higher resistance and resilience in mixed vs. pure stands, with very similar fluctuations in all the climate scenarios studied.
The complementary effects hypothesis links tree species diversity and stand productivity to explain the advantages of mixed stands over monocultures (Aarssen 1997;Loreau 2000;Steckel et al. 2019;Van de Peer et al. 2018). According to this hypothesis, mixed forests with varying leaf areas, canopy heights and growth rates can obtain and utilize resources more effectively, thus improving stand productivity (Tilman et al. 1997). Jucker et al. (2014) argue that Pinus sylvestris stands let enough light pass through the canopy to allow broadleaf plants to grow below. Similarly, morphological differences in the canopies of Pinus sylvestris and oak complement each other and increase light utilization throughout the forest (Pretzsch et al. 2015a, b). For the mixed coniferous forests analyzed in the present work, the pure-stand accumulation rate of Pinus sylvestris, Pinus pinaster, and Pinus nigra was very similar (Fig. 4), and their accumulation rate curves ran almost parallel to each other. Mixed stands of Pinus sylvestris with either of the other two conifer species also showed a relatively similar accumulation rate that decreased by nearly the same value (about 50%) by the end of the simulation, indicating overyielding of about 28% compared to pure Pinus sylvestris stands.
Of the combinations studied, our results indicated that Pinus sylvestris-Fagus sylvatica and Pinus sylvestris-Pinus nigra mixed stands had smaller differences with respect to pure stands of Pinus sylvestris, Fagus sylvatica and Pinus nigra. Previous studies have reported higher accumulation rates for Fagus sylvatica in mixed stands compared to pure stands (Del . This was expressed in terms of overyielding (Condés et al. 2013) and even transgressive overyielding (Preztsch and Schütze 2009) and attributed to species niche complementarity (Del . 2017). Our results for Pinus nigra were quite similar, although differences between pure and mixed stands were lower during all the periods. Specific studies for Pinus sylvestris-Pinus nigra interactions are lacking in the literature, but since P. nigra is a close relative of P. sylvestris (Gernandt et al. 2005), the similar behavior and niches of both species could explain the scant differences between pure and mixed stands. In a recent study, Aguirre et al. (2019) found that mixed stands of two Pinus species generally had a neutral or negative effect on each other, resulting in similar or lower productivity than in pure stands. Toïgo et al. (2015) also suggested that overyielding would occur in some cases and underyielding in others, depending on the degree of complementarity and competition between the two but also on abiotic conditions. They also observed that combining two conifers does usually not lead to significant overyielding, which is broadly in line with our simulation results.
In contrast, overyielding was seen in Quercus pyrenaica -Pinus sylvestris mixtures at the end of the study period. This result is consistent with previous studies of Pinus sylvestris -Pinus pinaster (Riof. 2017a, b) and Pinus sylvestris -Quercus pyrenaica mixtures (Muñoz-Gálvez et al. 2021). In conifer-conifer mixtures, vertical structure complementarity due to differences in growth velocity could lead to higher light-use efficiency, while in the conifer-broadleaf mixtures, complementarity in crown and root systems may imply higher belowground resource efficiency (Forrester et al. 2014;Forrester and Bauhus 2016). However, although the results for the two mixtures indicated overyielding, both were found to benefit stand regeneration in prior studies (Del . 2009;López-Marcos et al., 2020). This suggests a positive effect on long-term stand stability.
One of the main findings of the present work is that the mixed stands in the study showed a smaller decline over time. In the trend that emerged, mixed stands eventually outperformed or widened the gap in relation to pure stands. In the first stage, the CO 2 accumulation rate remained high and stable or increased slightly in all four mixtures (Fig. 4). This may have been related to initial forest conditions. (SDI was not at maximum, competition between trees was less intense, and forests were growing faster.) In our simulations, we established a 2% reduction for all trees when SDI reached a constant SDI max . Theoretically, this would allow for relatively high growth until the forest reached SDI max . Then, after SDI > SDI max , in pure stands, CO 2 growth would enter a multi-year cycle featuring a short sudden increase followed by a decrease, and then stability [regardless of climate, i.e., De Martonne aridity index (M)]. In mixed stands, each 2% reduction had the potential to change the species composition, driving the forest to the most stable species ratio. In this scenario, the growth rate in mixed stands would moderately exceed that of the stable period for pure stands. We did not reach this stage in the 100-year timeframe of the simulation. However, the curve for mixed forest CO 2 production rate presented a small, sudden fluctuation at the horizontal coordinate of the year 2020 due to the fact that M did not vary continuously in our simulations but is updated every 20 years (Fig. 4).
Our results may also suggest that the various combinations of tree species analyzed may respond differently to harsh living conditions in future. In the simulations, M decreased with SSP, so that higher M values were found in SSP 1 (most optimistic scenario) than in SSP 5 (least optimistic scenario). Similarly, the total CO2 stock for all pure and mixed stands decreased as SSP became more severe (decreasing M), and M began to act on growth and mortality. However, mixed-stand accumulation rates were significantly higher than those of their respective pure stands under the studied scenarios.
In our findings, overyielding in terms of CO 2 accumulation rates appeared in all four species mixtures. However, there was a clear difference between the results obtained for the mixture of Pinus sylvestris with the two broadleaf species. The CO 2 accumulation rates in the Pinus sylvestris 1 3 -Fagus mix were 30-37% higher than in a pure Pinus sylvestris stand. In contrast, Pinus sylvestris -Quercus mixes produced only 24% overyielding relative to pure Pinus sylvestris stands, which was the lowest result of the four groups. Although the accumulation rate of pure Fagus stands declined rapidly as the study progressed, its high initial value may suggest that it could rapidly occupy a portion of the mixed stands. This means that the mixed forest would not be dominated by either species and both would contribute to overproduction. The same would apply to the Pinus sylvestris -Quercus mixture. The simulations clearly indicate that Fagus sylvatica is well suited for mixing with Pinus sylvestris. Reality also bears this out, as the two species are complementary in the shape of their canopies and the spatial use of their root systems (Pretzsch et al. 2015a, b). Pinus sylvestris roots usually occupy only the upper 0-40 cm of soil, whereas Fagus sylvatica roots can reach depths of 40-80 cm and can help Pinus bring water from deeper to shallower layers. Fagus also brings minerals from the deeper layers to the surface through its own growth cycle, in the form of leaf litter, etc. This improves the humus layer (Pretzsch et al. 2015a, b) and upper soil minerals, optimizing the soil environment for Pinus. However, drought tolerance in Pinus-Fagus stands is a concern that may limit future viability. Martín-Gómez et al. (2017) have indicated that Quercus trees may have a significant survival advantage over Pinus trees under prolonged drought conditions. The decline of Pinus forests and their subsequent conversion to Quercus forests has been observed in many locations. The same was observed in our simulation results, where the CO 2 accumulation rate for Quercus declined only moderately between SSP1 and the more arid SSP5. The Pinus-Quercus mixture could therefore be much more drought tolerant than the Pinus-Fagus mixture.

Conclusions
Innovative simulations of the evolution of CO 2 sequestration capacity in Mediterranean pure and mixed stands were developed and presented in this work using the SIMANFOR platform and data from the Spanish National Forest Inventory. These simulations predicted the evolution of various pure and mixed forests under four climate scenarios. Higher CO 2 stock and accumulation rate values accompanied the most optimistic climate change scenario (SSP1) and lower values were detected in the most pessimistic scenario (SSP5).
Pure and mixed stands differed in terms of CO 2 stock and CO 2 accumulation rate over time. At the beginning of the simulation period, mixed stands generally presented lower CO 2 stock than their respective pure stands, but the differences were drastically reduced by the end of the simulation period in favor of the mixed forests, as reflected in the CO 2 accumulation rate.
Interest in conifer-broadleaf mixed stands has increased in recent years, and forestry policies can be designed to support species mixtures in the effort to create more resilient stands. Our results showed similar behavior in the conifer-conifer combinations, but the conifer-broadleaf mixtures exhibited opposite trends. While Pinus sylvestris -Fagus sylvatica had the lowest accumulation rate reduction, Pinus sylvestris-Quercus pyrenaica had the highest. Mixed forests in the study exhibited overyielding compared to pure forests, which highlights the importance of mixing species that will yield higher CO 2 values at the end of the simulation period.
Further research is needed to reduce vulnerability in Mediterranean forests through better understanding of how mixed stands evolve under different management options. Replicating simulations of this kind can inform the development of silvicultural guidelines and measures that will help forest ecosystems adapt to climate change.