Species interactions under climate change in mixed stands of Scots pine and pedunculate oak

Mixed-species forests have become widely studied in the recent years because of their potential to mitigate risks associated with climate change. However, their growth dynamics are often difficult to predict because species interactions vary with climatic and edaphic conditions, stand structure and forest management. We examined species interactions in mixtures of Scots pine (Pinus sylvestris) and pedunculate oak (Quercus robur) under climate change and for varying soil conditions in the Netherlands, over a period of 30 years. We parameterized, calibrated and validated the 3-PGmix model for mixing effects in Scots pine and oak mixtures and analysed these effects under climate change. 3-PGmix performed well for the variety of forest stands examined throughout the Netherlands. Furthermore, it was also able to reproduce mixing effects for each species in mixtures compared to monocultures for the growing conditions examined. Simulated climate change resulted in lower productivity of oak and higher productivity of Scots pine, compared to the current climate. This was observed for both monospecific stands and mixtures. The mixture of Scots pine and oak showed clear but limited overyielding (mixture yield greater than the mean of the monocultures), which was mainly attributed to oak. This was maintained under the most extreme climate scenario for 2050, implying that for oak, increased growth due to mixing with Scots pine was larger than the reduction in productivity under the future climate. On resourcelimited soils, Scots pine competitiveness was increased, and this was maintained under a warmer and drier climate. Our results suggest that projected changes in climate will influence species interactions and result in increased Scots pine productivity, notably on poor sandy soils, which are typical of the Netherlands.


Introduction
Climate change is expected to increasingly impact forest ecosystems worldwide. Adapting forests to climate change is therefore one of the major challenges for forest management. For this, it is necessary to understand the effects of climate change on forest systems, as well as to be able to predict how these effects might change over time (Keenan, 2015). It is generally assumed that higher temperatures can facilitate growth processes (Saxe et al., 2001), and higher growth rates for temperate and boreal forests have been reported in many studies (Boisvenue and Running, 2006;Kauppi et al., 2014;Pretzsch et al., 2014;Reyer et al., 2014;Bussotti et al., 2015;Pretzsch et al., 2018). However, prolonged periods of drought may hamper this growth, reduce vigor and resilience (Bréda et al., 2006;Hartmann, 2011;Pretzsch et al., 2013a;Buras et al., 2020) or even result in additional mortality (Allen et al., 2010;Anderegg et al., 2012;Rigling et al., 2013;Buras et al., 2018).
Mixtures are considered increasingly important in the face of climate change. Over the last few decades, the share of mixed forests in temperate regions has increased, and currently more than two-thirds of Europe's forests are dominated by two or more tree species (Forest Europe, 2015). Mixed forests provide more resilience and resistance to various biotic and abiotic climate-related disturbances, such as pathogens, herbivores, drought, wind and fire (Bauhus et al., 2017). Species interactions in mixtures may result in higher productivity (Vilà et al., 2013;Pretzsch et al., 2015a;Jactel et al., 2018) and forests with more tree species may also provide higher levels of multiple ecosystem services simultaneously (Gamfeldt et al., 2013). If the productivity of mixed-species stands exceeds the weighted mean productivity of the monocultures, this is called overyielding, or transgressive overyielding when it exceeds the productivity of the most productive monoculture (Hector et al., 2002).
The drivers of mixing effects have been discussed in many studies in terms of competition, facilitation and competitive reduction (Forrester and Bauhus, 2016). These species interactions in mixtures are often not static, and change in space or over time as resource availability changes, during stand development, or in response to climate change (Forrester, 2014). It is therefore crucial to understand under which conditions species interactions might change for a certain species combination, in the face of changing growing conditions due to global change. Since the impacts of climate change may differ between species and hence may also vary with species composition, mixing proportions and (a)biotic conditions, the direction of interspecific competition may change under climate change. For this reason, it is of the utmost importance for forest managers to understand species interactions and to possibly adapt management practices in these mixed stands, in the face of climate change.
Here we studied species dynamics of Scots pine (Pinus sylvestris L.) and pedunculate oak (Quercus robur L.) under different growing conditions and climate change scenarios for the Netherlands. These species dominate many Dutch forests: Scots pine dominates on about one third of the total forest share and oak is the most important broadleaved species (Schelhaas et al., 2014). Most forested land in the Netherlands is confined to poor sandy soils, and the mixture of pine and oak is one of the few species combinations that occurs on these soils (Den Ouden et al., 2010). Both species are of great importance for forestry in the Netherlands: Scots pine is valued for its timber that serves many purposes, ranging from construction wood to paper and pulp, and oak for its high quality timber and its role in nature conservation (Jansen et al., 2018a). Although Scots pine is considered only slightly more light-demanding than oak (Niinemets and Valladares, 2006), they differ in leaf phenology (evergreen vs deciduous) and in light-interception, which is higher for oak than for Scots pine (Bréda, 2003;Sonohat et al., 2004;Balandier et al., 2006;Forrester et al., 2018). Some studies on mixing effects for Scots pine with either pedunculate oak, sessile oak (Q. petraea) or a combination of both, report overyielding, mainly driven by oak. These results are based on permanent plots throughout the Netherlands (Lu et al., 2016), with increased overyielding on poor soils (Lu et al., 2018), and a modelling study in France (Perot and Picard, 2012). Transgressive overyielding for Scots pine and oak was reported in two recent studies that analysed triplets of a mixture and two corresponding monocultures; in Germany and Denmark (Steckel et al., 2019) and in one third of the stands throughout Europe along a productivity gradient (Pretzsch et al., 2020a). No overyielding effects, however, were observed for the mixture in French NFI plots, although non-significant higher productivity in mixture was seen in oak (Toïgo et al., 2015a;Toïgo et al., 2015b). These studies suggested that overyielding resulted from light-use complementarity, due to differences in phenology, shade-tolerance and crown architecture. Bello et al. (2019) showed that contrasting water-use strategies of oak and Scots pine in mixture reduced water constraints during a drought period, notably for oak.
Climate-sensitive models designed for mixed stands are useful tools for providing insight into interspecific interactions, and associated physiological responses, under different growing conditions and climate change scenarios. The process-based forest growth model 3-PG (Landsberg and Waring, 1997), is such a model. It has been used to study climate effects on growth, is widely used in many different regions, under many climatic and site conditions, and has been parameterized for many species (Landsberg and Sands, 2011). It has been recently transformed from a model for even-aged monocultures of evergreen species into a model that is also applicable for mixtures and deciduous species, 3-PG mix (Forrester and Tang, 2016).
The aim of this study was to analyse mixing effects in stands of Scots pine and pedunculate oak under a changing climate. The first objective was to parameterize, calibrate and validate the 3-PG mix model for growth and yield of Scots pine and oak under Dutch growing conditions. The second objective was to reproduce mixing effects in mixtures of Scots pine and oak. Next, we investigated if species interactions in those mixtures would be different under climate change and how this may vary with soil conditions. For the analysis, a time-horizon of 30 years was used, to match available plots, to exclude regeneration, and to stay close to the horizon for management decision making. We hypothesized that growth and yield of both species, as well as mixing effects on productivity, would increase in all climate scenarios, although less in the most extreme scenario. Furthermore, we hypothesized that, under climate change, potential shifts in dominance of one species over the other would be observed first in the resource-limited soils.

Description of the model
The process-based stand level model 3-PG was developed by Landsberg and Waring (1997), and adapted for mixed-species purposes by Forrester and Tang (2016), hereafter called 3-PG mix . The stand is divided into species or cohorts. These can interact within the canopy to influence light absorption. This results in vertical gradients of radiation, vapour pressure deficit and aerodynamic conductance, which influences species interactions in terms of water use and soil water availability. 3-PG mix has a monthly time step and consists of 5 sub-models, which will be described below. Details on the model and the adaptations to the mixed-species version 3-PG mix can be found in Forrester and Tang (2016).
The first sub-model predicts light absorption in the canopy and calculates gross primary production (GPP) based on maximum potential light-use efficiency of the species (α cx ). The light absorption depends on species-specific light extinction coefficients, leaf area indices and the vertical positioning of the crowns. Light-use efficiency α cx is reduced in response to limitations resulting from non-optimal temperature, frost, vapour pressure deficit (VPD), soil moisture, soil fertility, atmospheric CO 2 and stand age (Landsberg and Sands, 2011). NPP (Net Primary Production) is calculated following a fixed proportion of GPP: NPP/ GPP = 0.47 (Waring et al., 1998). The second sub-model allocates NPP to the different biomass compartments of each species: foliage, stem and roots (McMurtrie and Wolf, 1983). Partitioning to roots is influenced by soil fertility, VPD and soil moisture. Tree size influences partitioning between stem and foliage. Sub-model 3 adjusts the number of trees per hectare using the −3/2 self-thinning law in order to account for density-dependent mortality. The forth sub-model calculates the soil water balance and the fifth sub-model converts biomass to output variables such as species-specific mean diameter, height, basal area and wood volume.

Relative productivity
A wide range of concepts and methods is used to characterize the productivity of mixed stands . Here, we chose to express mixing effects in terms of relative productivity (RP), since it accounts for species proportion (Williams and McCarthy, 2001 (2) in which p 1 and p 2 represent the growth or yield in monocultures of species 1 or 2 respectively, p 1(2) the productivity of species 1 in mixture with 2, and m 1 or m 2 the mixing proportions in terms of basal area. Values of RP higher than 1 indicate overyielding, and values below 1 underyielding, as shown in Fig. 1. In this conceptual diagram, the mixing effect in terms of relative productivity (RP) for Net Primary Productivity (NPP) is demonstrated for the mixture of species 1 and 2.
In the left panel, relative productivity of both species is above 1, but overyielding declines with increasing availability of nutrient resources, while in the right panel species 1 is overyielding, species 2 is underyielding and the mean of both stands (RP total stand ) indicates overyielding.

Field data collection and study site
In order to parameterize, calibrate and validate the model with field data, we utilized the long-term growth and yield data that were collected and maintained by the Forest Ecology and Forest Management Group of Wageningen University (Den Ouden and Mohren, 2016). This database is derived from permanent field plots in the Netherlands, and is a collection of growth and yield data of even-aged monocultures and even-aged mixtures of several species and 2-species combinations. For this study, the monocultures of Scots pine and pedunculate oak were used, as well as the corresponding mixtures (Table 1). Monocultures were defined as plots where at least 90% of basal area of the plot belonged to one species. In the mixed stands the species mentioned first is always the main species and the second the admixed species. The 3-PG mix model uses several input variables at the stand level, among which are planting date, initial number of trees and biomass components. As the stands were thinned regularly and the model is sensitive to tree density, thinning data, such as stem numbers of the remaining stand on each date of measurement, were included in the inputs for the model.
The Netherlands has a mild maritime climate, caused by predominant southwest winds, with moderately warm summers (average temperature 17 °C), cool winters (on average 3.4 °C) and high humidity. The mean annual temperature in the reference period 1981-2010 was 10.1 °C. Annual precipitation, on average 851 mm per year, is generally evenly distributed throughout the year (KNMI, 2014).

Estimation of biomass components
The 3-PG mix model requires biomass input for 3 different components: stem (including branches and bark), leaves and roots. Therefore, the model was initialized with the estimated amount of biomass on the first measurement date of the field plots. In order to do so, we used the allometric equations developed by Forrester et al. (2017a) to estimate initial biomass. We calculated individual tree biomass components for the different compartments (Y) using species-specific equations that include effects of diameter at breast height (DBH, B in Eq. (3)) and both DBH and stand basal area (BA in Eq. (4)) in the following functional forms: From these individual tree biomass data, we estimated initial Fig. 1. Conceptual diagram of the mixing effects in terms of relative productivity (RP) for species 1, species 2 and total mixture in relation to a certain abiotic factor. Values above y = 1 indicate a higher relative productivity in mixture than in monoculture and those below y = 1 a lower relative productivity (Eqs. (1) and (2)). In the left panel (a) both species indicate overyielding, whereas in the right panel (b) species 2 has a relative productivity below 1, indicating underyielding. biomass for stem, roots and foliage (Ws, Wr and Wf in Mg per ha, respectively). To overcome the problem that for some of the oak plots initial measurements were carried out in winter time, foliage biomass for those plots was estimated in the first month after leaf formation in spring.

Meteorological input data
The 3-PG mix model requires monthly mean values of temperature (Tav, in °C), maximum temperature (Tmax, in °C), minimum temperature (Tmin, in °C), monthly precipitation (in mm month −1 ), the amount of frost days (days month −1 ) and monthly mean daily solar radiation (MJ m −2 day −1 ). For these variables, we used homogenized weather data from the main climate station De Bilt in the Netherlands (KNMI, 2017), which has the longest time series of Dutch weather data. We used this station, because nearly all plots were situated closer to this weather station than to any of the other long term stations. Moreover, in the Netherlands there are no mountains influencing local weather conditions, and the country is relatively small. The weather data were available from 1901 onwards; except for data on precipitation that could be obtained from 1906 onwards and on solar radiation that were available from 1957 onwards. For those plots that were first measured before 1957, extrapolation for the missing intervals was carried out, utilizing monthly averages.

Other initialization data
Apart from stand characteristics such as stem number and age, biomass estimates, and meteorological requirements, the 3-PG mix model requires initialization of soil water balance and fertility.

Soil data
Soil data from the ISRIC-World Soil Information SoilGrids250m database (http://www.isric.org/explore/soilgrids) were used to determine the plot soil texture at a resolution of 250 m. This database is a system for automated soil mapping based on spatial prediction methods using soil profile and environmental covariate data. Predictions for sand, silt and clay content at seven standard depths (0, 5, 15, 30, 60 100 and 200 cm) were used. We calculated average sand, silt and clay content for each of the depth intervals (0-5, 5-15, 30-60, 60-100 and 100-200 cm) by averaging the predicted values of the upper and lower boundary of the interval. Then this average value was multiplied with the thickness of the depth interval in order to obtain weighted means. This resulted in a weighted sand, silt and clay content for every plot (Lu et al., 2018). With this soil texture, it was possible to determine the required input parameters describing soil texture (c θ and n θ ) and the potential available soil water (ASW), following Landsberg and Waring (1997).

Soil fertility
Soil fertility is quantified in the model by the fertility-dependent growth modifier (Fertility Rating, FR). It takes values between 0 and 1 for low and high soil fertility, respectively. Growing conditions for a forest stand, including soil fertility, are reflected in the site index (SI). We assumed that the SI of a forest stand depends on climate, soil water availability and soil fertility. This can then be described by the following function (Forrester et al., 2017b): where the variability ε that is not explained by climate or soil water availability is assumed to be a measure of the soil fertility. In this function, Martonne is the aridity index according to De Martonne (1926):

Aridity index
annual precipitation (mm) mean annual temperature ( C) 10 (6) and ASW is the potential available soil water, which is the soil water holding capacity based on soil texture and depth. Since SI and soil fertility represent long-term average conditions, a long-term De Martonne index was calculated for De Bilt weather station. The dominant height at age 50 years was taken as a proxy for SI per plot, using yield tables and the OPTAB model (Jansen et al., 1996). In order to determine soil fertility, we fit equation (5) with the known values (SI, ASW and De Martonne index) and assumed that FR can be expressed as the remaining variability in the predicted SI, which was not explained by soil water characteristics or climate. Therefore, we assumed that FR is the observed SI plus the residual, and the residual is the observed SI minus the predicted SI (Eq. (5)). Finally this fertility data was normalized, by dividing the FR estimation per plot by the maximum FR estimation per species, to obtain an index from 0 to 1 for each plot.

Parameterization
The description of the parameter values and their sources are indicated in Table A1. For Scots pine, the parameter set used in this study was adapted to the Dutch circumstances from the 3-PG set for Scots pine in monocultures in Scotland (Xenakis et al., 2008) and the European 3-PG mix set for Scots pine-beech mixtures (Forrester et al., 2017b). Most of the 3-PG mix parameters for pedunculate oak were based on field data and published literature. Other parameters were default settings (Sands, 2010) and remaining parameters were fitted. The latter comprise the parameters that could not be estimated reliably from the literature and that vary between species, thereby preventing the use of default values. These parameters quantify foliage to stem partitioning (p 2 and p 20 ), litterfall rates (γ Fx and γ R ) and root partitioning (η Rx and η Rn ).

Allometric relations in 3-PG mix
The canopy structure is calculated using allometric relationships that predict mean tree height, crown length and crown diameter. Mean height (h, in m) is calculated as a function of mean stem diameter at breast height (B, in cm) and competition (C), as described by the following equation (Eq. (7)): in which a H , n HB and n HC are constants. The competition index (C) is calculated as the sum of the species-specific products of basal area and basic wood density using Eq. (8) (Forrester and Tang, 2016): where BA is the basal area (m 2 ha −1 ) and ρ is the basic wood density (g cm −3 ) of species i (Zanne et al., 2009). The competition index is then expressed per ha after dividing by plot area (m 2 ). The equations for height were fitted as mixed models to estimate the parameters, with tree and measurement record nested within plot as the random effect of the models. The input data for these models were species-specific data from both monocultures and mixtures, in order to capture all the variance in height for the specific species in the Netherlands.
For our study sites, detailed information on mean live-crown length (hL, in m) and mean maximum crown diameter (K, in m) was not available, so we used parameter values from Forrester et al. (2017b) for Scots pine, and estimated them based on the functions for hL and K used in Forrester et al. (2017b) with data from long-term Swiss plots for pedunculate oak (Forrester et al., 2019). We used a correction factor to account for the bias that occurs when back-transforming logarithmic data (Snowdon, 1991).
The parameters that relate stem biomass and diameter (a s and n s ) were taken from Forrester et al. (2017a), as well as the parameters to relate age and specific leaf area. Values for mean specific wood density were used from Forrester et al. (2017b) for Scots pine and from Zanne et al. (2009) for oak.

Quantum efficiency and biomass partitioning
The model predicts GPP from light absorption and the maximum canopy quantum efficiency (α Cx ), and then converts this GPP to NPP and derives volume growth. It is assumed that maximum volume growth rates of each species can be taken as a base to estimate the α Cx parameter. These maximum growth rates, presumably without limitations imposed by temperature, frost and other factors, were obtained from maximum current annual increments in the Dutch yield tables (Jansen et al., 2018b, Jansen et al., 2018c, on the best sites in the Netherlands. These growth rates were then converted to NPP using biomass expansion factors and wood density (Vande Walle et al., 2005), then to GPP and this was divided by the absorbed light (APAR) based on estimates in Forrester et al. (2018). The light extinction coefficient (K H ) for Scots pine was based on Forrester et al. (2017b) and derived from Bréda (2003) for oak.
Parameters for the modifiers on α Cx such as temperature limits (T min , T opt , T max ), maximum stand age and effects of vapour pressure deficit on stomatal conductance (K D ), were taken from Forrester et al. (2017b) for Scots pine. For oak, stomatal conductance (K D ) was set at the default value and temperature limits and maximum stand age were approximations based on San-Miguel-Ayanz et al. (2016). The parameters that determine leaf production and leaf fall for the deciduous species (leaf P and leaf L ) were set at the 1st of April and the 1st of November, based on local knowledge. The maximum proportion of rainfall interception by the canopy (I RX ) was calculated from Augusto et al. (2002) and Van Nevel (2015) for oak and were taken from Forrester et al. (2017b) for Scots pine. Leaf litterfall rate γ F1 was calculated for Scots pine (Schoettle and Fahey, 1994;Oleksyn et al., 2003), following Landsberg et al. (2003) and fitted for oak.

Model performance: calibration and validation
Model calibration was executed with data from monoculture plots. The model was initialized with stand characteristics such as trees per hectare, stand age and estimations of stem, foliage and root biomass at the first date of measurement; and then run until the last date of measurement for each plot. Model predictions were compared with field data at the last date of measurement per plot, and hereafter each parameter was changed, one at a time for several rounds, until predictions matched observed data as precisely as possible. Since the model was calibrated using the monocultures, a validation was carried out with the mixed plots. Hereafter, the final parameter set was obtained and overall model performance was tested.
To test if the model reproduced mixing effects in terms of relative productivity (RP) that matched the observed field plots, we used corresponding plots with similar soil, climate and other site-specific conditions that may affect tree growth. The corresponding plots cannot be considered true triplets, consisting of a mixture and two corresponding monocultures in close proximity that were established for the purpose of studying mixing effects, as these were not available in the dataset in sufficient numbers. Plots that were identified as a triplet were not more than 30 km apart, were in the same age category of young (up to 40 years), medium aged (up to 100 years) and old (above 100 years) stands, and in the same or similar Dutch soil class on scale 1:50000 (https://www.pdok.nl/viewer/). More details are provided in Table A2.
Comparisons of predicted output variables and stand characteristics were tested on the following criteria: relative average error (e%, Eq. (9)), relative mean absolute error (MAE%, Eq. (10)), mean square error (MSE, Eq. (11)) (Janssen and Heuberger, 1995;Vanclay and Skovsgaard, 1997) and the model efficiency (EF, Eq. (12)) (Loague and Green, 1991), where O i and P i are the observed and predicted values, respectively, and O and P are the means. For the first 3 criteria, values closer to 0 (perfect fit), indicate better model predictions. Positive and negative values of e% (Eq. (9)) indicate over-or underprediction, respectively. Model efficiency (Eq. (12)) can be less than 0, which means that the model prediction is poorer than one resulting from merely using the mean, and can have values up to 1, indicating a perfect correlation between predicted and observed values. All statistical analyses were performed in R 3.6.1 (R Core Team, 2019).

Projections in future climate
The 'KNMI '14 Climate change scenarios' for the Netherlands (KNMI, 2014) were used for growth projections under future climate change. These are two warm (W H and W L, plus 2 °C in 2050) and two moderate scenarios (G H and G L , plus 1 °C in 2050), in which the H subscript indicates a high change in atmospheric circulation patterns, resulting in wetter winters and drier summers. The L subscript indicates a relatively low variability in circulation patterns, leading to smaller changes in precipitation in both seasons. Input data for 3-PG mix for the 2050 time horizon were generated with the 'KNMI Time series transformation tool' (Bakker, 2015). This tool made it possible to apply the change of the mean values and the variability as prescribed per climate scenario (W H , W L, G H and G L ) to the given historical range of temperature, precipitation and global radiation of the reference climate . The climate data of the years 2000 until 2015 were used as current climate.

Scenario testing
The calibrated model was then utilized to analyze the influence of future climate change scenarios, soil conditions, thinning intensity and stand age on the mixing effects. Therefore the model was initialized with the mean data of the observed plots at the first measurement and then run for 30 years. Apart from the means of foliage, root and stem biomass, initial stocking and stand age, also three thinning scenarios were designed: no thinning, moderate thinning from below and heavy thinning from below. This was done by adapting every 5 years the residual stocking (trees ha −1 ) and the biomass pools after thinning, following the Dutch growth and yield tables (Jansen et al., 2018c;Jansen et al., 2018b). The model runs started in 2020 with 45 year-old-stands of 3 species compositions, two monocultures and a mixture, with 7 levels of soil fertility, 5 climate scenarios, 3 thinning intensities and 5 soil water holding capacities, that represented the range of conditions in the observed plots. This resulted in a factorial design of 3 × 7 × 5 × 3 × 5 = 1575 model runs for a period of 30 years.

Model performance
The 3-PG mix model produced accurate predictions of diameter, height and basal area (Fig. 2) and yields of total and stem biomass (Fig. 3) of sampled field data. Predictions were carried out for monocultures and mixed stands of Scots pine and pedunculate oak in the permanent field plots in the Netherlands at the final date of measurement. Predicted yields correlated very well with the observed measurements in the monocultures (Figs. 2 and 3, R2 > 0.95) and mixtures (R 2 > 0.9 for basal area and stem mass, R 2 > 0.98 for all other compartments).
The comparisons between predicted and observed stem and total biomass and the variables directly derived from them were slightly better for monocultures than for mixtures, since the criteria for model performance showed lower mean values for monocultures (mean e%= 0.16, MAE%=14.56, MSE = 390.6, N = 151) than for mixtures (e% =4.39, MAE%=18.52, MSE = 374.6, N = 36). Mean model efficiency was higher for monocultures than for mixtures (mono EF = 0.67, mix EF = 0.56) for those variables ( Table 2).
Predictions of mixing effects in terms of relative productivity (Eqs. (1) and (2)) of the biomass components and their derived variables were generally accurate as well (e%=2.73, MAE%=20.6, N = 16) with a mean model efficiency of 0.56 (Table 2). In contrast, model predictions were less precise for the smaller foliage and roots biomass compartments and are not shown here.

Projecting growth in climate scenarios
Model projections in different climate scenarios from 2020 until 2050, initialized with the mean data of the Dutch permanent field plots, showed that Scots pine had a slightly higher total biomass growth under the future climate than the current climate (Fig. 4). This effect was smaller in monocultures than in mixtures. Oak showed the opposite pattern; lower biomass growth under the future compared to the current climate.
Furthermore, Scots pine total biomass in mixture, when corrected for mixing proportions, was lower compared to pure stands. This effect was observed in both the current and future climate.

Projecting mixing effects in climate scenarios
Simulated mixing effects in terms of relative productivity (RP) on stem biomass and NPP (Net Primary Productivity) were highest for pedunculate oak and above 1 (Fig. 5). In other words, this means that oak mixed with Scots pine showed a higher relative biomass than when growing in monoculture. The overall mixing effect was also above 1, indicating that this Scots pine-oak mixture was slightly more productive than the mean of its corresponding monocultures.
Scots pine, however, was less productive in mixture, as was also demonstrated in Fig. 4. There was a positive effect of future climate on Scots pine productivity, but the negative effect of mixing with oak had a much larger impact on Scots pine growth. In analogy with this, these projections show no reductions in productivity under a future climate for oak in mixture, while the positive effect of growing with Scots pine had a greater influence on oak growth than the future climate.

Projecting mixing effects along different gradients and thinning scenarios
As available soil water decreased, mixing effects for both Scots pine and pedunculate oak tended to decrease. For oak, this was observed in all climatic scenarios including the current climate. Scots pine, however, showed a different pattern: in the driest soils mixing effects were higher under the future W H climate than under the current climate (Fig. 6a). On those dry soils, Scots pine showed slight overyielding in terms of NPP in the most extreme climatic scenario (Fig. 7a).
Along gradients of soil fertility, mixing effects for both species diverged; differences in the effects were very small in the poorest soils but increased as soils became richer. Meanwhile, the effect of the future climate on mixing effects was small but slightly increased in the richer soils (Figs. 6b and 7b).
As the stands aged, oak stem mass mixing effects increased, but NPP aging effects were not clearly observed (Figs. 6c and 7c). The effect of thinning, however, was more pronounced: more intense thinning scenarios increased the differences between both species, which were further augmented under future climate change (Figs. 6d and 7d).

Validation
The 3-PG mix model was able to predict growth of stands of Scots pine and pedunculate oak between the first and final date of The model reproduced mixing effects consistent with other studies. In many studies that analyze productivity of mixtures compared to monocultures, triplets of two monocultures and one mixture were utilized (Pretzsch et al., 2013b;Pretzsch et al., 2016;Del Río et al., 2017;Forrester et al., 2017b;Zeller et al., 2017;López-Marcos et al., 2018;Steckel et al., 2019;Pretzsch et al., 2020a;Russo et al., 2020). Here we used more distant, but comparable plots in terms of soil and age, which we classified as triplets. We only included the most comparable plots in our database, which reduced the number of plots for measuring mixing effects considerably. Although this approach may have caused less accurate predictions compared to using triplets that were established especially for the purpose of studying mixing effects, the model error for mixing effects on stand characteristics was low (2.73%) and correlation between predicted and observed mixing effects above 90%, mostly above 97%. Unlike for empirical analyses, where it is important to use plots with similar ages, the age is accounted for by 3-PG, which makes it less important to have similarly aged plots. Moreover, if the plots are of different ages, it may even give a better indication that the model performed well. It is important to note that in the scenario predictions, triplets with exactly the same site conditions were used as model input.
The accuracy of the predictions of the stand variables was comparable to the European triplets of Scots pine and beech that were simulated with 3-PG mix (Forrester et al., 2017b), even though the sampling periods and timings, and therefore also the model runs, in the Dutch stand data were variable (mean = 16.7 years, min = 4 years, max = 57 years) and not fixed for 12 years from 2002 to 2014, as in the Scots pine and beech study. For example: some of the Dutch plots were measured from 1948 until 2002, while other plots were sampled from 1984 until 1996. Furthermore, over 90% of the plots were regularly thinned, with thinning intervals of 3-7 years. We wanted to include this variability in stand, soil, management and climatic characteristics in the model calibration, to be able to do predictions on as many forest stands as possible with the same parameter set.
The size and accuracy of the growth and yield data set used in this study was important for at least three reasons. Firstly, calibrations using plots along broad environmental gradients are effective at constraining parameters associated with water and nutrient sensitivity to a similar degree as nutrient, drought and irrigation experiments (Thomas et al., 2017). Therefore, the wide-spread distribution and high numbers of Table 2 Statistical information that describes relationships between predicted and observed variables for mixtures and monocultures of Scots pine and pedunculate oak. The equations for e% (relative average error), MAE% (relative mean absolute error), MSE (mean square error) and the EF (model efficiency) are described in the methods section, slope is of the relationship forced through the origin, P-value for the test of whether the slope of the relationship is significantly different from 1 and the R 2 values.  plots used in this study ensure that these parameter sets can be applied to a broad range of climatic and edaphic conditions. Secondly, the accuracy, in terms of information about stand age and recording which trees died as well as those that were thinned, improves the accuracy of the parameters. Other process-based model calibration studies found that experimental plot networks similar to ours improved the accuracy of calibrations compared with less accurate inventory data from National Forest Inventories (Minunno et al., 2019). Thirdly, our plot network included very long time series. It was expected that foliage and roots compartments would be less accurately predicted (Paul et al., 2007;Forrester et al., 2017b), because these compartments are relatively small in comparison to stem biomass, they have a larger variability and higher turnover time (Forrester et al., 2017a), and the biomass in the "observed" data was based on allometric relationships instead of being measured directly in the field. Therefore, emphasis was put on stem biomass and derived variables.
There are three notable sources of uncertainty in the input data. First, solar radiation measurements before 1957 were absent and therefore monthly averages were used for this period. Although the model starts with light absorption from solar radiation, we believe this could have led to minor uncertainties only, since the majority of the measurements were done after that date. The second source of uncertainty is the FR estimation. SI was calculated with the OPTAB model, which uses the widely applied function of Chapman-Richards, in the Dutch yield tables (Jansen et al., 1996) and therefore gap-filling approaches had to be done for dominant height at age 50. This could have led to higher uncertainties in the stands that deviate more from this age. As our field plots of both species were on average 48.9 years old, we do not consider this of major importance. This SI was then used as input to estimate FR, which is challenging to determine and considered to be the weakest part of the model, although the process of determining it has improved thanks to contributions of many different studies, as described in the review by Gupta and Sharma (2019). Biased FR estimates could have possibly led to biased model predictions, however, the results of our calibration and validation process (Figs. 2 and 3) do not give any reason to consider this to be a major source of uncertainty. The third source of data uncertainty is the use of allometric relationships for observed biomass values. We used the equations from the database of Forrester et al. (2017a) instead of direct measurements. Therefore, we not only took predicted and observed biomass compartments into consideration in the calibration process, but also stem diameter, height and basal area. In the model, these variables were namely derived from biomass, but in the observed data, they were directly measured in the field, herewith providing good insight in the performance of the model.
Several uncertainty and sensitivity analyses exist for 3-PG (Esprey et al., 2004;Xenakis et al., 2008;Song et al., 2013;Augustynczik et al., 2017) and also for 3-PG mix (Forrester and Tang, 2016). Although some small variations may be expected in the sensitivity analysis of each parameter set, the above-mentioned analyses show similar results. In general, the model is most sensitive to some of the least uncertain parameters (Xenakis et al., 2008), among which were quantum yield efficiency, specific leaf area, foliage to stem biomass partitioning ratio, the ratio of net to gross primary production, canopy conductance and basic wood density (Esprey et al., 2004;Xenakis et al., 2008;Forrester and Tang, 2016). These parameters impact many of the outputs, especially foliage biomass (Forrester and Tang, 2016). Root turnover rate was found to be the least certain parameter (Xenakis et al., 2008), but it was among the parameters that had the smallest effect on the output in our study. As the more sensitive parameter values were either taken from literature, or empirically-based, this ensured a good parameter fit. It is possible that the literature information used to calculate some of the parameters did not cover the range of provenances included in the Dutch populations where the model was applied. This could be improved as more studies on those processes are carried out. However, given the results of the validation, we do not consider this to be highly important.

Species interactions in Scots pine-oak mixtures
We observed higher mixing effects in terms of relative productivity (RP) for pedunculate oak in all gradients of environmental drivers. This was maintained over stand age and in future climate growth predictions, even though oak showed a growth reduction in the most extreme W H climatic scenario, while Scots pine showed the opposite. This indicates that the effects of mixing were larger than the climatic effects, since RP was always above 1 for oak and almost always below 1 for Scots pine, regardless of Scots pine's higher productivity under the future climate. The total mixing effect along environmental gradients was generally larger than 1. This overyielding has been observed in many other studies on mixtures of evergreen-deciduous species, e.g. Scots pine-beech (Condés et al., 2013;Pretzsch et al., 2015b), Norway sprucebeech (Pretzsch et al., 2010;Toïgo et al., 2015a) and Douglas-fir-beech (Bartelink, 1998;Lu et al., 2016;Thurm and Pretzsch, 2016;Lu et al., 2018), depending on stand characteristics and site conditions. For Fig. 5. Simulated mixing effects in terms of relative productivity of stem biomass (a) and of NPP (b) for Scots pine, pedunculate oak and total mixture. Patterns are shown for the year 2050 after 30 years of model run in each of the 5 climatic scenarios. Cur = current climate, whereas G L , G H , W L and W H are the KNMI '14 scenarios for 2050, arranged in order from less to more extreme, respectively. Values above y = 1 indicate a higher relative productivity in mixture than in monoculture and those below y = 1 a lower relative productivity (Eqs. (1) and (2)). mixtures of Scots pine and oak as studied here, overyielding was reported throughout the literature; in a modelling approach in one site of France (Perot and Picard, 2012), in empirical studies on the same Dutch plots as this study (Lu et al., 2016;Lu et al., 2018), in triplets in Germany and Denmark (Steckel et al., 2019) and along a productivity gradient in Europe (Pretzsch et al., 2020a). However, in a study using French NFI data (Toïgo et al., 2015a), the total mixture did not overyield, but on the species level, higher growth rates in mixture than in monoculture were observed for oak. This was confirmed in another study on the same plots (Toïgo et al., 2015b). In the European triplet study, overyielding was observed for many, but not all, Scots pine and oak mixtures (Pretzsch et al., 2020a). It was suggested that overyielding was driven by complementary light-use, due to differences in leaf phenology and crown architecture (Lu et al., 2016;Steckel et al., 2019;Fig. 6. Simulated mixing effects in terms of relative productivity of stem biomass for Scots pine, pedunculate oak and total mixture along gradients of soil water holding capacity (a), soil fertility (b), stand age (c) and thinning scenario (d). Solid lines are model projections in current climate, while dotted lines are projections in the most extreme W H KNMI'14 2050 scenario. Values above y = 1 indicate a higher relative productivity in mixture than in monoculture and those below y = 1 a lower relative productivity (Eqs. (1) and (2)). Mixing effects are shown for the year 2050 in all graphs, except for graph c. Pretzsch et al., 2020a;Steckel et al., 2020) and complementary wateruse efficiency (Brown, 1992;Bello et al., 2019).
It is challenging to compare these findings in the literature, as productivity, and therefore overyielding, can be defined and measured in many different ways (Sheil and Bongers, 2020) and analyzed using different methods, varying from empirical to modelling approaches.
Apart from issues related to definition and methodology, another possible explanation for the variation within and between studies, is the variability in climatic and soil conditions (Lu et al., 2018;Steckel et al., 2019;Pretzsch et al., 2020a). Although in the present study, higher mixing effects on productivity for oak compared to Scots pine were consistent throughout all environmental gradients, and overyielding Fig. 7. Simulated mixing effects in terms of relative productivity of NPP for Scots pine, pedunculate oak and total mixture along gradients of soil water holding capacity (a), soil fertility (b), stand age (c) and thinning scenario (d). Solid lines are model projections in current climate, while dotted lines are projections in the most extreme W H KNMI'14 2050 scenario. Values above y = 1 indicate a higher relative productivity in mixture than in monoculture and those below y = 1 a lower relative productivity (eqs. (1) and (2)). Mixing effects are shown for the year 2050 in all graphs, except for graph c. occurred for most site conditions, we observed differences in mixing effects on the poorest and driest soils and under future climate change.
Our model predictions showed that interspecific competition from oak resulted in lower productivity of Scots pine in mixtures compared to monocultures. This is consistent with other studies in temperate forests, where Scots pine mixing effects were smaller than those of oak (Perot and Picard, 2012;Toïgo et al., 2015a;Lu et al., 2018;Steckel et al., 2019). In the European triplet study, on the other hand, with plots along a variety of site and environmental conditions, overyielding was balanced between both species: in some plots oak was the better competitor while in others it was Scots pine (Pretzsch et al., 2020a). Similar findings were reported in a study in England (Brown, 1992).
We hypothesize that light-related competition from oak caused the reduced the growth of Scots pine. Several studies have suggested that this mixture has the potential for light-use complementarity (Lu et al., 2016;Steckel et al., 2019;Pretzsch et al., 2020a). However, light interception by Scots pine is likely to be lower than that of oak (Bréda, 2003;Sonohat et al., 2004;Balandier et al., 2006;Forrester et al., 2018). This lower light interception, which was accounted for in the parameterization of 3-PG mix , may thus have led to a generally lower light availability for Scots pine in mixture, which possibly could have led to reduced growth compared to pure stands. In addition, this may explain why interspecific competition experienced by oak was lower than intraspecific competition, as was also seen in another modelling study (Perot and Picard, 2012). Moreover, as oak is slightly more shadetolerant than the light-demanding Scots pine (Niinemets and Valladares, 2006), it has the ability to grow in different vertical strata of the stand (Perot and Picard, 2012;Lu et al., 2018;Steckel et al., 2019;Pretzsch et al., 2020a). Compared to monocultures, oak may thus be able to grow slightly overtopped by Scots pine (Pretzsch et al., 2020a).
Recently it was observed that Scots pine and pedunculate oak may increase height growth due to competition (Forrester et al., 2017c;Del Río et al., 2019), thereby potentially improving their access to light. In the present study, this change in allometry was included in the 3PG mix parameter n HC , which describes competition in the stem height relationship (Eq. (7)). Although the observed allometric effect of competition on height growth was slightly higher for Scots pine than for oak (Table 2), the predicted height differences between Scots pine and oak were small; in both the scenario simulations (< 2 m) and the field data (2.22 m). Thus, severe competition for light in the crown layer may have suppressed Scots pine growth due to its lower light interception and, in analogy with this, favored oak growth.

Productivity in climate change scenarios
In our climate change predictions with average soil water availability levels, oak showed reduced biomass growth in a hotter and drier climate compared to the current climate (Fig. 4). However, slightly increased productivity under the future climate was observed in monocultures of Scots pine and those effects were more pronounced in mixtures with oak. As oak is considered to be the more drought-sensitive species in the mixture (Albert et al., 2015), it may be expected that in the dry sandy soils in the Netherlands, future climate change towards the most extreme scenario has a higher impact on oak than on Scots pine. Indications that these growth reductions in oak were indeed driven by drought effects in the extreme scenario, are demonstrated by the fact that in the W L climatic scenario, which is a warm scenario but without the strong summer droughts, a small growth increase was observed in oak compared to the current climate (Fig. A1). This is in line with observed data on forest productivity in response to climate change in areas that are not limited by water stress (Boisvenue and Running, 2006;Pretzsch et al., 2014).

Mixing effects in climate change scenarios
As water availability increased, higher productivity in mixture at the stand level was observed (Figs. 6a and 7a), which was also seen in triplet studies (Steckel et al., 2019;Pretzsch et al., 2020a) and this is in agreement with a large meta-analysis on mixtures (Jactel et al., 2018). On the species level, however, the pattern is more variable. Although it is known that growth reductions may occur in Scots pine trees under drought stress (Zang et al., 2011), we observed augmented pine competitiveness on soils with low water availability under the more severe climate. This is in line with a Spanish study on mixtures of Scots pine with beech, that showed that interspecific competition had a positive effect on Scots pine growth in dryer sites (Condés and del Río, 2015). In contrast, Pretzsch et al. (2020a) found that pine was superior to oak at higher levels of water supply and the authors suggested that this was related to better adaptation of pine to colder and wetter conditions. Another recent study found that the drought response was different for both species and suggested that mixing may reduce the drought susceptibility of Scots pine and oak (Steckel et al., 2020). This may be explained by the results from Bello et al. (2019), in which Quercus petraea and Scots pine were found to utilize deeper water sources in mixture than in monocultures, hereby changing water availability in mixture in severe summer drought.
Furthermore, we found that Scots pine was most competitive in the poor soils (Figs. 6b and 7b), in line with another study that used the same plots (Lu et al., 2018). In contrast, mixing effects for oak productivity increased with higher soil fertility levels, which may be explained by its need for a higher nutrient supply (Fürstenau et al., 2007). In addition, we may reason that mixing improved the light conditions for oak, as increased mixing effects on richer soils are likely to occur when species interactions improve light-absorption or light-use efficiency (Forrester, 2014). For Scots pine the reverse may have been the case, since higher mixing effects on productivity on poor soils are likely to occur when species interactions improve soil resource availability (Forrester, 2014). We thus can hypothesize that, on these poor soils, Scots pine may have profited from released competition from oak, which is in line with the resource-ratio hypothesis (Tilman, 1985).
In the most severe climate change scenario, higher nutrient levels resulted in higher growth for both species (Figs. 6b and 7b). However, this was observed in combination with average soil water availability and not maintained in the driest soils for oak (Fig. A2). It can therefore be assumed, that with a hotter and drier climate, competitive strength may shift from oak to Scots pine, and even towards Scots pine dominance under the poorest and driest conditions. Altered competitive dominance in mixtures during periods of drought stress was also observed in other species (Cavin et al., 2013;Rubio-Cuadrado et al., 2018). As most of the mixed stands of oak and Scots pine occur on the poorest and driest sites in the Netherlands, the sandy soils in the central to eastern part of the country, we may expect the first effects of climate change there.
We have to emphasize that although 3-PG mix is a process-based model that accounts for the main ecological factors, it is a stand-level model. This implies that effects of climate change on individual trees cannot be simulated. For example, the model does not take into account that size may modify drought sensitivity (Zang et al., 2011;Bennett et al., 2015), nor does it include cumulative effects of drought periods, such as the consecutive dry summers of 2018 and 2019, which may affect some trees within the stand more than others. In addition, it must be emphasized that mortality due to other climate-related events such as storms, fires or insect outbreaks were not included here, although these effects can have large impacts on forest vitality, such as the major outbreaks of bark-beetle attacks in Norway spruce in central Europe (Pretzsch et al., 2020b). Finally, it is important to realize that the effects demonstrated in this study are based on model projections for 30 years, which is a relatively short period for oak-Scots pine forests. Therefore it is plausible, that the effects of global warming on species interactions will be larger in the long run.

Consequences for forest management
At the stand level, clear but limited overyielding productivity in terms of stem mass and NPP was observed for most of the environmental conditions and this was maintained under climate change. This emphasizes the potential of this mixture under global warming. At the species level, however, future climate may cause shifts in competitive strengths, in particular on dry and poor soils. As our results suggest that severe competition for light in the crowns may suppress Scots pine productivity in this mixture on the more humid soils and under a future warmer climate, increasing vertical stratification by creating unevenaged stands with Scots pine in the upper canopy could be beneficial. At the same time, although a shift towards Scots pine dominance may be foreseen on poor soils, it is expected that the mixture will maintain its relevance for forest management.

Conclusions
In this study, a forest growth model recently adapted for mixtures was used as a tool to help disentangle the factors influencing species interactions in Scots pine and oak mixtures and to allow projections under global warming. We observed that this mixture showed clear but limited overyielding in all climatic scenarios. Although this overyielding was mostly attributed to oak and was maintained under the future climate, our projections showed decreased oak productivity in a warmer and drier climate, hereby likely reducing its competitive strength. This effect was most evident on the poor and dry soils. Our results strongly suggest that, under global warming, the reduced competitive strength of oak on those soils may lead to an alteration in species competition in favor of Scots pine. The first signs of this shift would be expected on the poor and dry sandy soils in the Netherlands.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.