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 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. 2009, 2017; Pardos 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 (Durrant et al. 2016; 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 (Durrant et al. 2016), 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 CO2 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 and Schü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 (Del Río et al. 2016; 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, 2018; Del Río et al. 2017; González-de-Andrés et al., 2017, 2018; Pretzsch 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 CO2 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 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 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 (Alberdi et al. 2016; 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.

Table 1 Mean values (standard error in brackets) of the main variables for individual and paired species
Full size table
Fig. 1
figure 1

SNFI2 plot distribution of the four mixtures in Spain. a Pinus sylvestris – Pinus nigra, b Pinus sylvestris – Pinus pinaster, c Pinus sylvestris – Fagus sylvatica, d Pinus sylvestris – Quercus pyrenaica

Full size image

Data were taken from the WorldClim2 (Fick and Hijmans 2017) CMIP6 downscaled future climate projections. Mean annual temperature (T) and total precipitation (P) values from five global climate models (BCC-CSM2-MR, CanESM5, CNRM-CM6-1, IPSL-CM6A-LR, MIROC6) under four Shared Socioeconomic Pathways (SSP1, SSP2, SSP3 and SSP5) at four different periods (2021–2040, 2041–2060, 2061–2080, 2081–2100) were calculated and derived for selected SNFI2 plots using R programming (R Core Team 2020) and QGis software (QGis.org 2021). Data from period 2000–2020 were obtained from the same source, corresponding to real data collected instead of projections.

CO2 stock and accumulation rate simulation in Spanish mixed forests

In this study, using the SIMANFOR platform (Bravo et al. 2012), CO2 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.

Table 2 Tree growth equations model for individual and paired species
Full size table

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.

Table 3 Tree height-diameter equations for individual and paired species
Full size table

We controlled for mortality based on the maximum stand carrying capacity (SDImax), 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 SDImax 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 SDImax ≥ SDI, mortality does not occur in a specific plot for a given species, when SDImax < 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, 2012). 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 CO2 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:

Fig. 2
figure 2

Overview of Shared Socioeconomic Pathways (SSPs) representing combinations of mitigation and adaptation challenges related to climate change [Figure extracted from O’Neill et al. (2017)]

Full size image

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 energy-intensive, 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 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.

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 SDImax 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.

Results

Forest CO2 stock simulation for Pinus sylvestris mixtures during the 2000–2100 period

Simulated CO2 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 CO2 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 CO2 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 CO2 stock in the mixed stands compared to Quercus pyrenaica pure stands (9.7%) and lower CO2 stock compared to Pinus sylvestris pure stands (47.5%). Differences with the latter reduced in successive timespans to the point that CO2 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.

Fig. 3
figure 3

Total CO2 stock (Mg · ha.−1) simulation for the 2000–2100 period in the most realistic SSP scenario (SSP2)

Full size image

Aboveground and belowground CO2 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 CO2 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 CO2 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 CO2 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 CO2 accumulation rate simulation for the 2000–2100 period

Based on simulated CO2 stock values, accumulation rate in terms of megagrams of CO2 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 CO2 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 CO2 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 CO2 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 CO2 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.

Fig. 4
figure 4

Mean CO2 accumulation rates (Mg · ha−1·year.−1) for the different studied mixtures and tree fractions for the 2000–2100 period based on the intermediate scenario (SSP2)

Full size image

The CO2 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% CO2 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). Mixed-stand 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) CO2 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 CO2 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 CO2 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 CO2 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 CO2 stock and accumulation rates in mixed vs pure stands.

Our findings also indicated higher CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 Río et al. 2016). 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 out-performed or widened the gap in relation to pure stands. In the first stage, the CO2 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 SDImax. Theoretically, this would allow for relatively high growth until the forest reached SDImax. Then, after SDI > SDImax, in pure stands, CO2 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 CO2 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 CO2 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 CO2 accumulation rates in the Pinus sylvestris – 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 CO2 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 CO2 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 CO2 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 CO2 stock and CO2 accumulation rate over time. At the beginning of the simulation period, mixed stands generally presented lower CO2 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 CO2 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 CO2 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.