Modelling the influence of thinning intensity and frequency on the future provision of ecosystem services in Mediterranean mountain pine forests

Thinnings are silvicultural operations that repetitively reduce tree density to improve the vigour of the remaining trees and the economic benefits of the stands. Thinning can also enhance the provision of various ecosystem services throughout the rotation period. In this study, we use a modelling approach to evaluate stand dynamics and the provision of ecosystem services (mushroom production, wood for timber, carbon storage, blue water, and habitat for biodiversity) in Mediterranean mountain pine forests. We simulated thirteen thinning regimes—defined by intensity and frequency—under two different climate change scenarios. We implemented the thinning regimes using SORTIE-ND, an individual-based model of forest dynamics, and then we used models developed for the study area to calculate the provision of services. We used as a case study Pinus sylvestris and Pinus nigra stands, and we evaluated the effect of the thinning regime, climate change, and forest type over 100 years. Our simulations suggest that the production of timber and carbon decreased with increasing intensity and shortening frequency of thinnings, while the provision of mushrooms and blue water generally increased under those conditions. Large timber was better supported by a thinning regime with heavy intensity and shorter frequencies, which also favoured the rapid presence of large dead trees (dbh > 30 cm) in the stand. We found synergies between the production of timber and carbon, while the provision of mushrooms and particularly blue water were in a trade-off relationship with these services. Our simulations show that climate change will lower the production of ecosystem services under the most severe climate predictions and alter the effect of different regimes on the provision of services. We conclude that our modelling approach is a useful and efficient tool for answering questions that would otherwise require long-term studies, and that it can provide useful information to guide management efforts to adapt forest management to the challenges of climate change.


Introduction
Although traditional forest management has focused mainly on timber production, Mediterranean forests have always been characterized by their multifunctionality. The provision of ecosystem services has been gaining importance in forest decision-making and planning (Wolfslehner et al. 2019); thus, services such as carbon sequestration and fixation, nonwood forest products, biodiversity, and recreation are now included in the Criteria and Indicators for Sustainable Forest Management (MCPFE, 2003).
Integrating all the services that forests provide into forest management is a challenge, as different trade-offs and synergies between services often exist and must be considered (Duncker et al. 2012;Marchi et al. 2018;Morán-Ordóñez et al. 2020;Zanchi et al. 2014). Moreover, part of the ecosystem services can bring an economic use-value, either direct (e.g., timber, hunting, mushroom production), or indirect (e.g., watershed protection, soil conservation, carbon sequestration), while others bring non-use value (e.g., conservation of biodiversity, habitats, and species) (Croitoru 2007; Communicated by Thomas Knoke.

3
MedECC 2020). Stakeholders' preferences should thus be integrated into the decision-making process, as some of the services might be of more interest than others (Diaz-Balteiro et al. 2017).
Climate change is another challenge that affects all ecosystems and requires changes in the human way of acting (IPCC 2014(IPCC , 2021. The Mediterranean region is identified as one of the most vulnerable climate change hotspots, where increases in temperature and decreases in precipitations will lead to an increase in the frequency and severity of droughts and aridity (MedECC 2020). The main regional models predict significant increases in average annual temperature, with rises that by the end of the century could range between 3 and 8ºC depending on the emissions scenario considered (IPCC, 2021). The consequences of climate change on forest health, dynamics and functioning are already noticeable and are expected to worsen in the future, jeopardizing the provision of ecosystem services (Allen et al., 2010;Lindner et al. 2010;Morán-Ordóñez et al. 2021).
Adaptive forest management can contribute to adapting forests to future climate conditions. Some recent works have evaluated how different forest management options can affect the provision of ecosystem services (Ameztegui et al. 2017;Díaz-Yáñez et al. 2021;Mina et al. 2017;Morán-Ordóñez et al. 2021). Among the strategies proposed to foster climate change adaptation, thinning is considered a particularly useful tool (Vilà-Cabrera et al. 2018). Thinnings are applied to improve the productive and protective functions of the stand by reducing the competition among trees. They promote the growth of the remaining trees Torras & Saura 2008) and lessen the negative effects of water scarcity by reducing stand density (Ameztegui et al. 2017;Fernández-de-Uña et al. 2015). The effect of various thinning regimes on stand structure has been profusely studied in the field , as has the effect of forest management on the provision of one or more ecosystem services. However, the effects of different complete thinning regimes on forest dynamics and the simultaneous provision of multiple services have hardly been investigated. It is, however, necessary to better understand synergies and trade-offs between services and to guide decision making.
In this study, we adopt a modelling approach to explore the effect of thinning intensity and frequency on forest dynamics and the provision of ecosystem services (mushrooms, timber, carbon, blue water-i.e., fresh surface and groundwater-and habitat for biodiversity) in Mediterranean pinewoods. We implemented the thinning regimes using the SORTIE-ND model, which provided us with annual information about stand dynamics (composition and structure). We then integrated the provision of ecosystem services using models available for our study area. We explored the effects of thinning intensity and frequency over 100 years, and under two different climate change projections -RCP 4.5 andRPC 8.5 (IPCC 2014, 2021). Considering that forest practitioners cannot wait to observe the long-term effects of their management decisions, our study proposes a solid tool to forward evaluate the impacts of the decisions made today on the state of forests in the future. The integration of service provision in planning under different scenarios allowed us to generate valuable information that can contribute to climate change mitigation and adaptation in Mediterranean forests.

Site characteristics
We chose for our study two of the most common forest species in the Pyrenees and Pre-Pyrenees region in Catalonia (north-eastern Spain). Pinus sylvestris covers 10% of the total forested area. It is mainly distributed in the Pyrenees and Pre-Pyrenees Mountains, where it finds the optimum conditions at elevations between 900 and 1,600 m a.s.l. (Piqué et al. 2014). On the other hand, Pinus nigra covers 6% of the total forested area in Catalonia and is distributed in the Pre-Pyrenees Mountains and some areas in central Catalonia (Piqué et al. 2014), at elevations between 400 and 1000 m (Piqué et al. 2014). Pinus nigra is a light-demanding species; however, Pinus sylvestris has even stronger light demands (Beltrán Barba et al. 2012;Piqué Nicolau et al. 2011).
For this study, we created theoretical stands located at the average topographic and climatic conditions for each species in Catalonia, assuming this would correspond to their ecological optimum. The Pinus sylvestris stand used for this study is situated at an elevation of 1200 m in a mesic site, where the mean annual temperature is 7.6 °C and the annual precipitations are on average 860 mm (Table 1). The Pinus nigra stand is situated at an elevation of 700 m in a mesic site, where the mean annual temperature is 12.7 °C and the annual precipitations are on average 650 mm (Table 1).

Climate change scenarios
We simulated two scenarios of future climate according to the projections of Representative Concentration Pathways (RCP) 4.5 and 8.5, defined in the IPCC's Fifth Assessment Report (2014). RCP 4.5 represents a climate scenario based on the assumption that carbon emissions will be reduced by the end of this century, while pathway 8.5 assumes business as usual conditions where carbon dioxide emissions will be sustained until the end of the century. Under the RCP 4.5 projections, the global temperature will increase by 2.3-2.9 °C until 2100, whereas the RCP 8.5 scenario predicts an increase in the global temperature of 4.1-4.8 °C (IPCC 2014). Climate data were based on the predictions developed by the EU-CORDEX project (available at Earth System Grid Federation; http:// esgf. llnl. gov/) for the period 2020-2100. Data were obtained according to the CNRM-CERFACS-CNRM-CM5 global model that was later regionalized at 11 km resolution. To acquire future climate conditions for each forest stand, we downscaled and biascorrected the regional projections to 1 km resolution using the R package "meteoland" . meteoland provides utilities to estimate daily weather variables at any position over complex terrains through spatial interpolation of daily weather records from meteorological stations and statistical correction of meteorological data series from climate models. A more detailed introduction to the package functionality can be found

Initial stand structure and thinning regimes
The initial characteristics of the stands were extracted from ORGEST, the regional guidelines for the management of Pinus sylvestris (Piqué Nicolau et al. 2011) and Pinus nigra (Beltrán Barba et al. 2012) in Catalonia. The stands were initiated from the age of 20 years, when the first thinning is proposed according to the ORGEST models. The number of trees per hectare (N) at the age of 20 years was 2500 for Pinus sylvestris stand (PiS) and 2100 for Pinus nigra stand (PiN), while the mean diameter of the stand was 10.7 cm for PiS and 9.6 cm for PiN. The diameter distribution of the initial stand was defined using a Weibull distribution function. The parameters for PiS were N = 2,500 trees/ha, Dg = 11 cm, and standard deviation (sd) = 3, while for PiN, N = 2,100 trees/ha, Dg = 10 cm, and sd = 3 (Table A.1).
For each species, we tested several thinning regimes, which consisted of different combinations of intensities and frequencies. The intensities considered were light (extraction of 15% of the stand basal area), moderate (25%), and heavy (35%). In all cases, thinnings were applied from below, i.e., we removed the trees from smallest to largest until the target basal area was reached. This type of thinning is typically applied in Spain to shade intolerant species such as pines, to reduce competition by removing mostly dominated trees, while also reducing vulnerability to wildfires . Each of these intensities was applied at different time intervals (frequencies)-every 10 years, 15 years, 20 years, and 25 years until the age of 50 to 70 years depending on the frequency. Thus, 13 different thinning regimes were defined for each species, including a non-intervention (control) regime. The details of each thinning regime can be found in Table A.2.

Modelling the thinning regimes and forest structure dynamics
Forest dynamics was implemented using SORTIE-ND version 7.05.07 (http:// www. sortie-nd. org). SORTIE-ND is an individual-based and spatially explicit model that simulates forest dynamics, as well as the effects of forest management and other disturbances (Canham et al. 2004). SORTIE-ND combines a series of empirical and mechanistic processes to determine the demography of individual trees in a stand, including their growth, reproduction, and mortality. Stand dynamics result from the combination of processes affecting each tree, which are species-specific and spatially explicit, as each tree has a specific location within the plot (Cristal et al. 2019).
The parameters needed to run the simulations were obtained using likelihood methods and data from the Spanish Forest Inventory from previous works (Ameztegui et al. 2015;Gómez-Aparicio et al. 2011) and have been successfully validated and used in previous studies in the area (Ameztegui et al. 2017;Morán-Ordóñez et al. 2020). SOR-TIE-ND is particularly appropriate for simulating the effects of thinnings on forest dynamics, as tree growth depends on tree size, climate, and the amount of competition exerted by their neighbours (Canham et al. 2004), which is estimated from the Neighbourhood Competition Index (NCI). NCI increases with the size of the neighbours and decreases with their distance from the target tree so that any change in the neighbourhood of a tree has an immediate effect on its growth (Eq. 1).
The calculation of NCI sums over the j = 1…n neighbours of each target tree i, where DBH j is the diameter at breast height of neighbour j; dist ij is the distance from target tree i to neighbour j, and α and β are parameters to estimate by the model. Details on the exact formulation of tree growth and the parameter values used to run the simulations can be found in Appendix B in Supplementary material. (1) The thinning regimes were implemented using the "harvest" behaviour in SORTIE-ND. This module allows the user to define any harvest operation by setting up (i) the year of application; (ii) the species it will be applied to; (iii) the cut type (i.e., partial cut or clear cut); (iv) the intensity of the harvest (in terms of percentage of stem density or basal area to remove); and (v) the type of harvest (from below or from above). Several harvest operations can be applied in a single simulation, and their combination defines the thinning regime that will be implemented.
All the details needed for one simulation are set through a parameter file which contains details about the plot (number of timesteps, years per timestep, location and climate), initial stand characteristics, and the behaviours or modules used during the simulation. An example of a parameter file used in thus study can be found in the Supplementary Materials. Overall, we created 130 parameter files for each species (13 regimes × 2 RCP = 26 × 5 repetitions = 130). The repetitions were considered to include the stochasticity of processes such as initial tree position-which is set randomly by the model-and mortality.

Modelling ecosystem services
SORTIE-ND produced one output file for each combination of species, RCP, thinning regime, repetition, and timestep, which contained observations for individual trees regarding the type of tree (alive, dead, harvested), and its position, DBH, and height. Using this information we modelled, for each output file, the provision of several ecosystem services: (i) mushroom production, (ii) timber production, (iii) carbon storage, (iv) provision of habitat, and (v) provision of blue water. We estimated mushroom production using the model developed by Bonet et al. (2010), where the annual mushroom production depends on the basal area of the stand, topography (slope, elevation, and aspect), and the total rainfall during autumn-August, September, October, and November. We calculated timber volume using species-specific allometric equations as implemented in the AllometrApp of the Catalan Forestry Lab (https:// labor atori fores tal. creaf. uab. cat/ allom etrapp). The overall amount of timber provided by each regime is therefore the volume of harvested timber after each thinning plus the final volume at the end of the rotation (100 years). We calculated the amount of carbon stored from the total amount of biomass-both above and belowground-existing at each timestep. The equations for each biomass compartment-which includes bark, wood, crown, and below-ground biomass-were obtained from the allometric equations included in Allometrapp (Ameztegui et al. 2022). Then, we calculated the amount of carbon from the amount of biomass by applying the conversion coefficient of 0.509 (Montero et al. 2005). We expressed border using the abundance of large dead trees-with DBH bigger than 30 cm (Cordonnier et al. 2014). Finally, we estimated the amount of blue water produced by each regime by coupling the outputs of SORTIE-ND with the water balance model included in the medfate R package (De Cáceres et al. 2015), as in Ameztegui et al. (2017). medfate performs daily updates of soil water content as a function of the stand structure and daily weather (radiation, temperature, and precipitation). Soil water balance is the difference between processes determining water input (precipitation) and outputs (canopy interception, tree transpiration, bare soil evaporation, surface runoff and deep drainage). We calculated annual estimates of blue water as the sum of runoff and deep drainage, and relative blue water was the ratio between blue water and precipitation for a given year. A more detailed description of the equations used for each ecosystem service in this study is presented in Appendix C.
Finally, we calculated, for each regime, the accumulated provision of each service over the rotation period (2020-2100). We assessed the effects of thinning intensity and frequency-as well as their interaction-climate scenario, and tree species on service provision using MANOVA, after testing the assumptions of normality and homogeneity of variance using the Kolmogorov-Smirnov and Bartlett tests, respectively. We assessed the significance of each explanatory variable on the combination of all the ecosystem services, and on each of them separately, via Wilks' Lambda, and we used partial eta squared as a measure of effect size associated with each main effect and interaction in the MANOVA model. We evaluated trade-offs and synergies of ecosystem service provision by computing pairwise Pearson's correlation coefficient between overall projections of ecosystem services provision. All data import, management, analyses, and visualization were performed in R version 4.1.0.

Effects of thinning regimes on forest dynamics
For both species, thinning regimes with a heavy intensity achieved the largest diameters at the end of the rotation period (Fig. 1a), while shorter frequencies resulted in smaller final diameters. The removal of 35% of the basal area every 10 years resulted in diameters twice as large as those of the non-intervention regime. However, the application of heavy and frequent thinnings comes with a cost, since it implied the extraction of more than 100% of the current growth. Thus, when these regimes were applied, the basal area (BA) maintained low throughout the rotation period, reaching values below 10 m 2 /ha in Pinus nigra stands and below 20 m 2 /ha in Pinus sylvestris stands (Fig. 1b). However, when applying low-intensity thinnings, the basal area quickly recovered after the interventions, reaching 30-40 m 2 /ha at the end of the rotation (Fig. 1b). The highest final BA values were obtained for the control plot, although regimes with long frequency and low intensity reached similar values. Final stand diameter and basal area were always lower under the climate predictions of RCP 8.5, although the differences between climate scenarios were greater for Pinus nigra (15-20% for d g; up to 30% for BA) than for Pinus sylvestris (5% for d g and 10% for BA).

Mushroom production
When a non-intervention regime was simulated, the combined effect of basal area and precipitation was clearly observed in our results. Thus, the highest production of mushrooms occurred when the stand reached around 30 years. After this, it decreased mostly due to the high Fig. 1 Effect of thinning frequency (every 10, 15, 20, and 25 years) and thinning intensity (removal of 0, 15, 25, and 35% of basal area) on the dynamics of tree mean diameter a; and basal area b during a rotation period under two climate change scenarios (RCP 4.5 and RCP 8.5) in one Pinus nigra (PiN) and one Pinus sylvestris (PiS) stand 1 3 basal area accumulated, but with slight increases in production for years with rainy autumns. Conversely, when different thinning regimes were simulated, the production of mushrooms was favoured after each thinning intervention, regardless the precipitation of that particular year (Fig. 2). In general, the highest mushroom production resulted when the heaviest intensity was simulated (extraction of 35% of basal area), but there was a strong interaction with thinning frequency: if frequent and heavy thinnings are applied, the basal area drops below the optimal (around 20 m 2 /ha) and results in declines in mushroom production. Interestingly, the differences between thinning regimes decreased in the most severe climatic scenario, indicating a major role of climate as compared to management. For example, for P. nigra under the RCP 8.5 scenario, there were almost no differences in mushroom productivity across thinning intensities and frequencies, and the total provision was very similar for all regimes and to the control simulation (Fig. 2). Climate change seemed to have a lower effect on mushroom provision in the Pinus sylvestris stand, although the production decreased anyway in the last 30 years of the simulation, when autumn precipitations are predicted to decrease under RCP 8.5 (Fig. 2).

Timber production
The simulated total timber production increased with frequency and, particularly, with decreasing intensity of thinnings, although the differences between intensities were reduced for longer frequencies. The highest total timber productions were obtained for regimes with low intensity (15%), which were more productive than the control plots on all occasions (Fig. 3). In the case of Pinus sylvestris, the regimes with moderate intensity were slightly more productive than the non-intervention regime, and with long frequencies, they were as productive as those of light intensity regimes. Pinus nigra, however, showed a lower response capacity against moderate or heavy thinnings. In the case of Pinus sylvestris, the most productive regime consisted of light intensities (15%) applied at intervals of 15 to 20 years. However, the difference with the control plots was not as great as for Pinus nigra, and the thinnings of moderate intensity already resulted in less timber volume than in the non-intervention regime (Fig. 3). Productivity was always lower in the most severe climate scenario, although the differences between scenarios were more notable for Pinus nigra (30%) than for Pinus sylvestris (10%).

Carbon storage
The highest total amount of carbon stored was reported for the control treatment-around 100 t/ha (Fig. 4). However, although the carbon stored was reduced after each intervention, the amount of total carbon stored by the end of the simulation was similar for the regimes with light thinnings and for the non-intervention plot. On the other hand, simulating strong and frequent treatments caused an abrupt decrease in stored carbon, which could not be recovered through post-thinning growth. In the most extreme cases, the carbon accumulated at the end was even lower than that available before the first thinning. The negative effect of climate change on the amount of carbon stored increased as the stand approached the end of the century, with a stronger effect in the case of Pinus nigra than for Pinus sylvestris.

Provision of blue water
Despite the interannual variations associated with meteorological fluctuations, the simulations predicted a marked effect of both thinning frequency and intensity on the provision of blue water (Fig. 5). In a non-intervention regime, less than 40% of the incoming precipitation would be released into aquifers, and there was a decreasing trend with time. This trend was reverted with management, but only when more than 15% of the basal area was removed. In the case of the heaviest thinning regimes-10-year frequencies and 35% intensity-the production of blue water sustainedly increased until the time of the last intervention, reaching 80% of incoming precipitation for Pinus nigra, and > 60% for Pinus sylvestris. Thinning frequency also had a role in blue water production, although lower than intensity, and there was a clear interaction between both factors. Indeed, the differences between intensities almost disappeared for the longest frequency, particularly in the case of Pinus sylvestris. There were clear differences in water use between both stands, with Pinus sylvestris releasing less blue water than Pinus nigra. However, there were no marked differences between both climatic scenarios, probably because we were using the relative amount of blue water as referred to incoming precipitation (Fig. 5).

Presence of large dead trees
The highest number of large dead trees (DBH > 30 cm) was obtained in the high intensity treatments (Fig. 6). However, this response was not linear, since more wood extraction often meant less mortality of the remaining stand, thus reducing the number of standing dead trees. In general, the availability of large deadwood decreased with the severity of the climate scenario. The reason was not a decrease in mortality, but rather a decrease in growth so that fewer trees reached a diameter greater than 30 cm.

Trade-offs and synergies of total ecosystem services production along a rotation
All the independent variables considered in this study (intensity and frequency of thinnings, climate scenario, and forest type) significantly and strongly affected the total production of ecosystem services throughout a rotation period (p-value < 0.05), and the size effect of all the explanatory variables was large ( 2 p > 0.5; Table D1). We also observed a strong effect of interactions between thinning frequency and intensity, and between each of them and the climatic scenario. In general, the total production of timber and carbon decreased with increasing intensity and thinning intensity (removal of 0, 15, 25, and 35% of basal area) on the number of large dead trees with DBH higher than 30 cm dur-ing a rotation period under two climate change scenarios (RCP 4.5 and RCP 8.5) in one Pinus nigra and one Pinus sylvestris stand and decreasing time between interventions, while the production of blue water and mushrooms-only for P. sylvestris-increased under those conditions (Fig. 7). In the case of Pinus nigra, the highest simulated total provision of mushrooms corresponded to high thinning intensity but moderate frequencies (15 years or longer). These intermediate treatments also seemed to favour the creation of large snags. Nevertheless, all the thinning regimes favoured the presence of large dead trees in comparison with a control plot.
We found a strong positive relationship (synergies) between the provision of timber and carbon, with Pearson's Fig. 7 Total production of ecosystem services (mushrooms, wood, carbon, blue water, and large dead trees) over a rotation period for a Pinus nigra and Pinus sylvestris stand under climate change projections RCP 4.5 and 8.5 r > 0.9 in all scenarios ( Figure D2 from Appendix D). On the other hand, we found a strong negative relationship (trade-off) between these two services (the provision of timber and carbon sequestration) and the provision of mushrooms and blue water, while the presence of large dead trees showed more moderate correlation with other services.

Discussion
Using a modelling approach, we could infer how thinning intensity and frequency affect the provision of various ecosystem services under two contrasted climatic scenarios. This kind of analysis can be useful in adapting forest management and planning to climate change, where the expected provision of ecosystem services can be integrated into the decision-making process. Overall, we observed that a thinning regime with light intensities and long frequencies provided an overall higher production of carbon and timber, and ensured more large dead trees to serve as habitat for species. However, a thinning regime with heavy and frequent harvests favoured a higher diameter growth and a higher provision of blue water, although all the other ecosystem services were at lower rates. The production of mushrooms was in turn favoured by intense, but more spaced-out harvests.

Effect of thinning regimes on forest dynamics
One of the traditional objectives of thinning is to change forest structure, so that the productive potential of the stand is concentrated on a smaller number of trees, but with a larger diameter, and therefore a higher market value. Our model correctly reproduced the positive effect of thinning intensity on diameter growth, even when thinning from below is applied, as proven in previous observational studies (e.g., del Río et al. 2008;Marchi et al. 2018;Montero et al. 2001;Ruiz-Peinado et al. 2017;Varmola & Salminen 2004). On the other hand, the effect of thinning intensity on basal area or volume is not so clear. Ameztegui et al. (2017) and del Río et al. (2008) showed that low-intensity thinnings did not decrease significantly the final basal area in comparison with a non-intervention plot. Yet, excessively heavy thinnings may lead to a loss in forest stock when the growth of the remaining trees cannot compensate for the loss of the extracted ones, as we observed here. As for the frequency, it is well known that the effects of thinnings are transient and disappear around 20 years after the intervention (Ameztegui et al. 2017). Our results show that the final productivity of the stands only matches that of undisturbed plots when frequencies are long. Low final densities can be convenient when the objective is the production of good quality, large-size timber (Montero et al. 2001). However, it should be noted that in our study, the thinning intensity was kept constant during the whole simulated period, while Montero et al. (2001) recommend reducing intensity with stand age, particularly in the last thinning.

Effect of thinning regimes and climate on ecosystem services provision
Overall, thinning regime proved much more important than the climatic scenario in determining forest dynamics and ecosystem services. The main effect of climate change represented by our models was to lower the growth of the trees, as previously observed in the area, especially for Pinus sylvestris (Fernández-de-Uña et al. 2015;Gracia et al. 2002). The lower growth rate automatically lead to carbon and volume losses, as reported before (Alvarez et al. 2016), but also to a decrease in most services. The decline in service provision under RCP 8.5 was particularly important for the Pinus nigra stand, which we placed in a more xeric site, in accordance with its distribution in the study area. However, meridional Pinus sylvestris may also suffer, since this species is known to have problems at the edge of its distribution (Benito Garzón et al. 2008;Fernández-de-Uña et al. 2015;Morán-Ordóñez et al. 2021). Interestingly, climate showed a strong interaction with the thinning regime for most services, suggesting that the choice of the best treatment should be made according to the ecosystem service to be promoted, and considering the most plausible climatic scenarios in each study area.
The combined effect of basal area and autumn precipitations for mushroom production was more deeply studied in empirical studies (e.g., Collado et al. 2018;Sánchez-González et al. 2019). Our modelling approach under different thinning regimes allowed us to see the effect of these factors on a longer term. The models we used predict a non-linear response of mushroom productivity to changes in basal area. Accordingly, we found that the lowest production of mushrooms was in a non-intervention regime, where the high basal area accumulation prevented mushroom formation. These results agree with Bonet et al. (2012), who found that Lactarius, mushrooms highly appreciated in the Pyrenees, dropped its productivity by 11% in the absence of management. In this scenario, mushroom production was mostly coupled with climate, and showed a sharp decrease by the end of the simulation period, as predicted by recent studies in the area (Morera et al. 2022). Mushroom production was favoured by the application of heavy thinnings, but only if enough time was elapsed between thinnings. de-Miguel et al. (2014) found that the yield increases with thinning intensity, since a decrease in canopy cover mean more water available in the soil. However, an excessive decrease may foster evaporation rates, negatively affecting fungal production (Bonet et al. 2012). We also observed that the differences in mushroom productivity between the two stands simulated were greater in the most severe climatic scenario-Pinus sylvestris produced three to four times more the mushrooms of Pinus nigra, as predicted elsewhere (de-Miguel et al. 2014;Morán-Ordóñez et al. 2020). However, the differences between regimes almost disappeared under RCP 8.5. This would indicate that, in the case of severe climate change, the effects of management are overridden, and the main determinant of mushroom production is climate, in particular autumn precipitation (Karavani et al. 2018;Morera et al. 2021Morera et al. , 2022. This was particularly evident for the Pinus nigra stand, located in a drier environment and therefore more dependent on variations in autumn rainfall. Total timber production decreased with thinning frequency and, above all, with the intensity of the interventions, although the differences between intensities were reduced for longer frequencies, suggesting an interaction between both variables. In general, the highest total timber volume was obtained for weak intensities (15%). Predictions for Pinus sylvestris showed greater adaptation to strong thinnings, as shown by the fact that moderately intense regimes resulted equally or more productive than control plots, while for Pinus nigra moderate thinnings always produced lower wood volume than unthinned ones. The effects of thinning intensity on wood production are widely known (e.g., del Río et al. 2008Zanchi et al. 2014), while the effect of frequency has not been as well studied. Varmola & Salminen (2004) showed the importance of early thinning in optimizing timber production, especially of large dimensions. Although thinnings decrease the biomass at the stand level, they increase diameter growth ), thus increasing the economic value of the forest (Baumgras et al.,1989). Despite total yields can be similar, the larger tree size we obtained in the simulations with more intensive treatments could represent higher economic returns, since larger trees usually imply more profitable destinations and lower unit extraction costs (Montero et al. 2001). On the other hand, faster growth can lead to lower wood density, which could have repercussions on wood quality (Jaakkola et al. 2005;Peltola et al. 2007;Russo et al. 2019).
Many studies have shown that one of the direct effects of forest management is the reduction of on-site total carbon stock compared to a non-intervention regime (e.g., Duncker et al. 2012;Ruiz-Peinado et al. 2013Zanchi et al. 2014). When the off-site carbon is considered, the carbon stocked by managed forests is however usually higher (Bravo-Oviedo et al. 2015). Ruiz-Peinado et al. (2016) found that the total carbon stock-including in-situ carbon and ex-situ carbon from removed trees-was not significantly different between an unthinned plot and a plot with light thinning regime. Yet in our study, when a thinning regime with light intensity and long frequencies was simulated, the level of in-situ carbon at the end of the simulated period reached the level of a non-intervention forest, while the total carbon stock was highest in the thinning regimes with light intensity. The differences between treatments would have been greater if we had considered the ex-situ carbon fixed by the wood products resulting from thinnings.
The positive effects of thinning on the provision of blue water are largely known and constitute the basis of the socalled eco-hydrological forest management (del Campo et al. 2017). According to our simulations, thinnings have the potential to revert the negative trend in blue water provision expected with climate change, but only when heavy and frequent thinnings are applied. More interestingly, the positive effect of heavy thinnings disappeared if more than 20 years elapse between thinnings, particularly for Pinus sylvestris under severe climate change, which matches previous studies on the duration of thinning effects (Ameztegui et al. 2017). This can be related to the higher adaptation of Pinus sylvestris to heavy thinnings that we mentioned above, but also to changes in the ratio between blue and green water when water becomes scarce (De Cáceres et al. 2021).
There is a lack of studies showing the effect of thinning on habitat provision in Europe and especially in the Mediterranean region, but more suitable habitat is typically found in unmanaged forests (Torras & Saura 2008), or when management practices that require less intensive interventions are applied (Duncker et al. 2012;Tomao et al. 2020). Our model included the inbuilt self-thinning behaviour for each species, and consequently low-intensity thinning regimes or control resulted in a higher presence of standing dead trees. However, the slow growth in the control plots led to very few-if any-large dead trees. Large-diameter trees are especially important for biodiversity as they host species such as bryophytes, lichens, fungi, and saproxylic beetles, and provide high-quality nesting habitat for many bird species (Paillet et al. 2010). In this regard, our simulations showed that intermediate thinnings can produce the highest amounts of large size snags (Blaser et al. 2013), which mainly appeared in the second half of the rotation period, as the simulated stands mature and approached the end of the century.

Trade-offs and synergies among multiple ecosystem services
When forest management is oriented towards multiple uses, not all ecosystem services can be equally maximized. In the regimes we have defined, timber and carbon production were positively correlated. Synergies between these services have also been shown by Schwaiger et al. (2019), who found that current carbon sequestration is negatively correlated with timber production, while total carbon sequestration is positively correlated with it. Although we found a weak correlation of provision of habitat (expressed as large dead trees) with the other services, most of the studies have shown trade-offs between wood production and biodiversity (e.g., Duncker et al. 2012;Schwaiger et al. 2019;Sing et al. 2018), especially for high-intensity treatments. The provision of blue water and mushrooms were in general negatively correlated with the other services. This result is especially relevant in water-oriented silviculture, as it makes it difficult to integrate blue water provision in multiple-use management decisions. Similar results were obtained by Ameztegui et al. (2017) and Morán-Ordóñez et al. (2020) for the Mediterranean region, while Duncker et al. (2012) found that forest management had a low influence on water quantity under continental conditions. Considering the trade-offs and synergies between ecosystem services, forest management at the landscape-scale should be designed considering the demands of different stakeholders in each location, or with a diversity of approaches to maximize the overall provision of services (Schwenk et al. 2012;Sing et al. 2018).

Management implications
Pine forests in the region have multiple uses, such as wood and non-wood forest products, pasture production, fire prevention, hydrological protection, biodiversity conservation, and recreation. However, management recommendations are still mostly based on the obtention of certain tree sizes that ensure good profitability for timber. In this manner, for an even-aged Pinus sylvestris forest on a highquality site, where the goal is producing trees with diameters of ~ 50 cm, Piqué Nicolau et al. (2011) recommend a thinning regime from below with moderate to heavy intensities at every 10-20 years. For an even-aged Pinus nigra forest, from a high-quality site, where the goal is producing trees with diameters of ~ 40 cm, Beltrán Barba et al. (2012) recommend a thinning regime from below with moderate intensities every 10-20 years. In our study, we have shown that a thinning regime with heavy intensity every 10 or 15 years can lead to the diameter targets proposed by Piqué Nicolau et al. (2011) and Beltrán Barba et al. (2012), but this would automatically affect the provision of the other ecosystem services. Thus, we do agree with Montero et al. (2001) that an intervention regime that maintains a low density is justified only when the ultimate goal is large size timber and the stand is in a highly productive site. Other management recommendations by the same authors, which consider the protective function of these forests, advocate for less intense and frequent thinnings, which according to our results would be more compatible with a varied service provision. Although the role of multifunctional forest management often seeks to maximize the joint production of services, local aspects need to be considered, and the global optimization does not need to be the best management option. In fact, our study can complement current management recommendations by allowing us to quantify and integrate into the decision-making the provision of several services under different climatic scenarios. This can allow managers take the final decision with all the available information, and depending on the local context and the specific demands of society.

Limitations of the modelling approach
In this study SORTIE-ND resulted an efficient tool to simulate forest dynamics for pine forests parameterized in Catalonia. However, to be used for other regions and species it needs to be parametrized accordingly, which is a timeconsuming and challenging endeavour (Evans & Moustakas 2016). In addition, SORTIE-ND does not provide direct estimations on the provision of services, and other models need to be used. In this exercise we employ a wide range of service provision models, from purely empirical equations to process-based models. Depending on the complexity and the factors involved in such models, they will be able to detect the interactions between management, climate and local aspects, and the results will be varying in reliability. In addition, although we tried to include as many services as possible, the list is short when compared to the many services provided by forests: erosion, soil protection, cultural values, recreation, etc. Therefore, despite some ecosystem services are hardly quantifiable or modellable, future works should try to gather information on as many services as possible, so they can be integrated into forest management planning.
Finally, another limitation of our study is that we did not include disturbances of any kind in the simulations. However, we know that the most visible effects of climate change will be through an increase in the frequency and the intensity of disturbances (IPCC, 2021). Future research should explicitly explore the effects of different disturbance regimes associated to climate change scenarios.

Conclusions
Our study shows that in a forest managed towards multiple uses, the long-term planning of thinning operations can have a profound impact on reaching the goals for the forest. We found, using a modelling approach, that the intensity and frequency of thinnings significantly influences the total production of ecosystem services over a rotation period. Our simulations have shown synergies between the production of timber and on-site carbon, which were favoured by less intense and frequent thinning regimes. On the other hand, blue water and mushroom provision were in a trade-off with the rest of the services and were favoured by more intense thinning regimes.
As regards our modelling approach, we acknowledge that SORTIE-ND, although it comes in a package with some limitations, was able to capture the effects of the thinning regime on forest structure and dynamics, matching previous results from observational studies in the field. Therefore, we strongly believe that our modelling approach is a useful and efficient tool for answering questions that would otherwise require long-term studies, and that it can provide useful information to guide management efforts to adapt forest management to the challenges of climate change. In any case, the ultimate choice corresponds to forest managers or planners, who should adapt their decisions based on the goals, the possibilities that each forest can offer, and the local market or population's needs and demands.

Acknowledgements
The authors want to thank Lluís Coll for helpful discussions during the conceptualization of this work, and the MSc European Forestry programme (Erasmus Mundus Joint Masters scholarships) for their support to Diana Simon during her stay at the University of Lleida.
Author contribution AA and DCS conceived the research and designed the experimental setup. DCS performed the simulations and the statistical tests and prepared the figures. Both authors wrote and reviewed the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This project was funded by the Spanish State Research Agency through the UMBRACLIM project (PID2019-111781RB-I00).

Competing interests
The authors declare no competing interests.

Conflict of interest
The authors have no competing interests as defined by Springer, or other interests that might be perceived to influence the results and/or discussion reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.