Uncertainties related to climate change and forest management with implications on climate regulation in Finland

Uncertainties related to climate change and forest management with implications on climate regulation in Finland Jari Vauhkonen ⇑ , Tuula Packalen Natural Resources Institute Finland (Luke), Bioeconomy and Environment unit, Yliopistokatu 6, FI-80100 Joensuu, Finland article info Article history: Received 6 October 2017 Received in revised form 8 December 2017 Accepted 19 February 2018 Available online 3 March 2018 abstract Forests play an important role in one of the most important ecosystem services, climate regulation. In order to mitigate climate change, various international agreements aim at decreasing emissions through Land-Use, Land-Use Change and Forestry (LULUCF) activities. In a legislative proposal by the European Union, emissions from forests are accounted for in relation to an estimate of average emissions for a range of years in the past. However, different forest structures, management activities, growth variations and impacts of changing climate may result in considerably different future emissions. We assessed the magnitude of potential uncertainties due to changing climate and forest management to the projections of carbon stocked in aboveand belowground forest biomass in Finland until 2050. We used an areabased matrix model, which was developed to incorporate climate-induced tree growth as a timeinhomogeneous Markov chain. The potential amounts of both the carbon stored and extracted varied considerably depending on the level and allocation of future harvests. If realized, climateor management-induced growth improvements could increase the carbon stocks by up to one third in the end of the simulated period. Projections based solely on business-as-usual transitions and harvests could therefore lead to inefficient decisions regarding future carbon stocks and harvesting possibilities. Ó2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/). 1. Introduction Climate regulation is one of the most important ecosystem services (Costanza et al., 2017), for which forests and forestry play an important role. On one hand, carbon accumulates through growth of trees into forest growing stock. On the other hand, Land-Use, Land-Use Change and Forestry (LULUCF) activities impact carbon stocks. First, sustainable land-use and forest management can conserve or increase forest carbon stocks. Second, harvested woodbased products and their bi-products can replace fossil-based products, materials and energy. The complexity of trade-offs between trees left growing for climate regulation service or those harvested for provisioning services – and consequently for climate regulation as diverse products – place challenges to decision making regarding LULUCF activities and their regulation. The issues related to LULUCF are reflected by regulation measures in various international agreements under the United Nations Framework Convention on Climate Change (UNFCCC). In line with the Paris Agreement (UNFCCC, 2017), the European Commission (EC) has presented a legislative proposal (EC, 2016) to set a binding commitment for each member state to ensure that accounted emissions from land use are entirely compensated by an equivalent removal of CO 2 from the atmosphere through action in the sector, known as ‘‘no debit rule”. In the proposed rules, emissions from forests are accounted for in relation to a so called national forest reference level. The forest reference level is an estimate of the average annual net emissions or removals resulting from managed forest land within the territory of a member state. The average values are calculated for a range of years in the past, referred to as a reference period, which is a political decision that applies for all countries. In practice, the past forest growth, harvests and, consequently, net emissions may vary a lot between years within a country and the variation pattern over time is not the same in all countries. First, there is variation of sinks due to the growth variation of trees (Mäkinen et al., 2002). Second, there is interannual variation of emissions (Andersson et al., 2007), due to the market fluctuation, for example. In addition, climatic variation affects the tree growth and the impacts of changing climate are assumed to vary between regions depending on their ecological conditions (Charru et al., 2017) and management (Henttonen et al., 2017). Furthermore, global market turbulence accelerates fluctuation in roundwood market in Europe (Packalen et al., 2017) and the impacts of global https://doi.org/10.1016/j.ecoser.2018.02.011 2212-0416/Ó2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding author. E-mail address: [email protected] (J. Vauhkonen). Ecosystem Services 33 (2018) 213–224 Contents lists available at ScienceDirect Ecosystem Services journal homepage: www.elsevier.com/locate/ecoser changes are expected to realize at different phase in different countries depending on their socio-economic situation. Consequently, the use of the same reference period for all countries may penalize or benefit a country depending on how well the future forest resource and market upturn and downturn coincide with each other compared to the national reference value. Forest resource projections have been used to model future biomass supplies at the European and national levels. There is a large variation in data and models used for the national projections tailored for local conditions and information needs (Barreiro et al., 2016). For example, sustainability and the provisioning of other ecosystem services than biomass-based products are taken into account as a series of ecological, economic and technical constraints, which limit the availability or accessibility of forests for wood supply (Alberdi et al., 2016). The European studies are often summed up from national projections (Schelhaas et al., 2017) carried out using a pan-European data set and a generic model that cannot fully account for differences between countries. Consequently, the results from the studies summed up to the European level often show a large variation, usually claimed to be a result of inherent uncertainties related to data (e.g., Rettenmaier et al., 2010; Bentsen and Felby, 2012) and biomass estimation models (Neumann et al., 2016). Even if context-dependent changes in future land-use, forest management, and climate obviously have implications that propagate the projections as uncertainties, studies addressing these aspects are missing (Barreiro et al., 2016; Schelhaas et al., 2017). According to the legislative proposal of EC (2016), the future carbon pools for estimating the emissions should be projected assuming a ‘‘continuation of current forest management practice and intensity” to make the emission accounting comparable between other sectors and member states, but also enable accounting for country-specific forestry dynamics. Based on this principle, Grassi and Pilli (2017) described a simulation framework, which is (a) parameterized by the prevailing forest age-class structure, increments, and business-as-usual harvesting practices and intensity; and (b) used to project the carbon pools after the reference period, assuming that harvests are continued in a similar magnitude as in the reference period, but relative to the development of biomass available for wood supply (for details, see especially Box 1 in Grassi and Pilli, 2017). However, because of the multiple factors causing variation to the future scenarios, as reviewed above, the member states would most likely benefit from the assessment of future forest carbon sink uncertainties when negotiating on the national forest reference level. Bayesian inference techniques such as Markov chain models (e.g., Nabuurs et al., 2000; Thürig and Schelhaas, 2006; Eriksson et al., 2007; Verkerk et al., 2011) may be applicable for quantifying uncertainties (cf., Smith and Marshall, 2008), due to the potential to flexibly vary the assumptions related to future scenarios (see also Vauhkonen and Packalen, 2017). The aim of this study is to test a Markov chain model for assessing the degree of uncertainties in the projections of carbon stocked in aboveand belowground forest biomass in changing climate and in the context of LULUCF regulation in Europe. The main objective is to develop methodology to incorporate tree growth variation and impacts of changing climate into an area-based, Markov chain model developed for projecting different management scenarios (Vauhkonen and Packalen, 2017). The secondary objective is to apply the methodology to quantify the uncertainties related to forest carbon, and consequently, to the selection of national forest reference level in the LULUCF regulation for Finland. Our analyses fundamentally cover two ecosystem services: roundwood harvests as a provisioning service and the related effects on the carbon extracted and stored in the remaining growing stock as a regulation service. 2. Material and methods 2.1. Overview As the Markov chain model, we used v. 2.0. of the European Forestry Dynamics Model (EFDM), which is implemented in the R statistical modeling environment (R Core Team, 2016) and can be downloaded from the EFDM project repository (FISE, 2017) as open source under the European Union Public License (EUPL). The EFDM is an area-based matrix model, in which the matrices represent forest areas classified according to ecological and socio-economic factors. The EFDM simulates the development of the forest area distribution as a product of its initial state, proportions of areas expected to be managed according to different silvicultural practices, and the corresponding transition probabilities. The transition probabilities are conditioned on the activities, which can both differ between factors such as site type, species, owner, and other factors either affecting the forest dynamics or needed for reporting. As elaborated by Sirkiä (2012) and Packalen et al. (2014), there is a transition matrix per factor combination and per activity. The initial state and activity and transition probability matrices can be derived through a simple classification and aggregation routine from National Forest Inventory (NFI) plot data, if the future development is assumed to follow that realized in the past. Growth models or simulators modifying the pairwise observations can be used to derive transition probabilities under changing climate. The EFDM was parameterized for the current climate using transition and activity probabilities derived from permanent NFI plots as described in detail in the open-access article by Vauhkonen and Packalen (2017). The framework was extended to include effects of climate change and convert the outputs to carbon. The analyses carried out here aim at quantifying the degree of uncertainty occurring, when decisions are made according to the transition probabilities observed in the past, but changes to these transitions occur due to the climate or management improving the growth. The general framework and especially these changes are described below, but regarding details of the parameters and their effects to the output, the reader is referred to the paper by Vauhkonen and Packalen (2017). The simulations were carried out in 5-year periods, which correspond to the measurement interval in the permanent plot data. Altogether eight periods were simulated, i.e., the last year of simulations is around 2050, depending on the initial measurement year. It was assumed that the land-use restrictions determined the silvicultural system applied. The development of forests without restrictions was simulated according to even-aged management and age and volume as the axes of the matrices. Forests with restrictions on wood supply were simulated according to an uneven-aged management, where final fellings were replaced with thinnings from above and the simulations were based on stem number and volume matrices. Only the natural processes were simulated for forests not available for wood supply. The derivation of the initial data and transition probabilities corresponding to observations made from the permanent NFI plots are described in Sections 2.2 and 2.3.1. For the uncertainty assessment, expected effects of climate and adapting the transition probabilities due to these changes were modeled using ‘‘a transition probability drifter”, as described in Section 2.3.2. Three different harvesting targets and two alternative allocations of the harvests were applied, as described in Section 2.4. The harvest decisions were based on roundwood volume (measured in m 3 ), but translated to carbon (tonnes) using output coefficients as described in Section 2.5. 214 J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 2.2. Initial state space The sample plot data from the 11th Finnish National Forest Inventory (NFI11; altogether 51,827 forest plots measured in 2009–2013) were used to the estimate the initial distribution of forest in an area of altogether 21.28 mill ha on productive and poorly productive forest land in Finland. The growing stock was 2,234 mill m 3 , which is approximately 95% of the entire growing stock in Finland, and only areas with low importance for forestry were excluded. For the analyses, the forests were classified to matrices with axes corresponding to either age and volume or stem number and volume under even-aged or uneven-aged management, respectively. To produce an adequate amount of observations for the estimation of the transition probabilities, the class limits for both the volume and stem number were derived as the values of the 10th, 20th, ..., 90th and 95th quantiles of the pairwise observations made from the permanent NFI plots. The age classes were defined as 0, 5, 10, .. ., 120, 120+ years, the class interval of five years corresponding to both the measurements and the simulation step used in the analyses. The class limits for the continuous variables are presented as Appendix A. These matrices were derived separately applying the following, static land-use classes: (i) known land-use restrictions: forests available, forests with restrictions on availability, and forests not available for wood supply; (ii) forest ownership: private, public + other; (iii) site fertility: altogether, five categories corresponding to four taxation classes traditionally used in Finland + fifth class including all poorly productive forest land; (iv) dominant species: pine, spruce, deciduous trees. 2.3. Transition probabilities 2.3.1. Current climate The transition probabilities corresponding to the current climate were derived using pairwise observations from permanent plots of NFI11 (altogether, 11,987 plots), which were measured approximately five years earlier in the previous inventory (NFI10). Positive differences in the total volumes on plots with no treatments based on data that could be matched with certainty between the two subsequent inventories were recorded as the pairwise data. The estimated transitions therefore included only growth and not potential reductions due to calamities or natural disturbances, for example. The transitions due to management activities were based on simulations of their expected development. The forests affected by final fellings were forced to transit to the beginning of the even-aged rotation. A thinning simulator was implemented to derive pairwise observations before and after the treatments. The thinnings took place in the beginning of each simulation period and the growth of the thinned forests in that period was simulated applying the transition probabilities of forests not managed. 2.3.2. Adapting the transition probabilities according to the expected climate change We used the following workflow, called ‘‘transition probability drifter”, to include the effects of climate-induced forest growth in the Markov chain model. In our drifter, a growth trend with stochastic variation was modeled in four steps, the first two of which are related to predicting increased CO 2 and temperature under climate change scenarios and the latter two steps to using these values for predicting the resulting growth increment. 1. CO 2 . Greenhouse gas (GHG) emissions resulting from climate change were expected to develop according to the Representative Concentration Pathway (RCP) scenarios. We used three out of the four scenarios, ignoring RCP6.0, as its effects with respect to forest growth were practically similar to those of RCP4.5 based on our modeling approach. In the scenarios RCP2.6, RCP4.5, and RCP8.5, the level of ambient CO 2 is expected to rise from the current level of 350 ppm to 443, 487, and 541 ppm, respectively, by 2050 (Meinshausen et al., 2011). 2. Temperature. The Finnish Meteorological Institute has projected the annual mean surface air temperatures in Finland to increase by 2040–2069 depending on the level of GHG emissions realized in the different RCPs. The expected changes in the annual mean temperatures, relative to 1981–2010, are expressed as normal distributions with parameters for the different RCPs (Ruosteenoja et al., 2016). The value representing the change in temperature by 2050 was obtained as a random value from the RCP-specific distributions. The temperature corresponding to 10%, 50%, and 90% values of the distributions were computed for the growth trend modeling. 3. Growth trend. Matala et al. (2005, 2006) used predictions obtained from a physiological growth model, FinnFor, to describe the impacts of elevated temperature and CO 2 on tree growth. The models presented give a ratio (or Relative Scenario Effect, RSE v ) of the volume growth under the changing climate to that under the current climate. Separate models, controlled by stand density, competition, site fertility and current temperature sum, are presented for Scots pine, Norway spruce and silver birch. To implement the models in conventional growth simulations, Matala et al. (2005, 2006) also express RSE v as a shape effect SE, which is a ratio of the height and volume growth, and present regression models for SE. The models of Matala et al. (2005, 2006) are available only for mineral soils; however, Nuutinen et al. (2006) modeled SE also for peatlands. We used the CO 2 and temperature values obtained according to the description above in the models to derive RSE V and used this ratio as the growth trend for our analyses. We computed RSE v for both mineral and peatlands using the models of Matala et al. (2006), but adjusted the RSE v values for peatlands according to the ratio of SE for mineral soil (Matala et al., 2006) to SE for peatlands (Nuutinen et al., 2006). 4. Stochastic annual variation in tree growth. In addition to the growth trend, we wanted to include the variations in tree growth as observed in the past growth series (cf., Henttonen et al., 2017). We used the autoregressive moving-average models fit to a trend-cleared growth-index series observed from Finland from 1890 to 1988 (Pasanen, 1998), which take into account the positive autocorrelation between successive years and the positive cross-correlation between the growths of different tree species. To obtain series of growth trends with stochastic variation, the previous four steps were repeated 1000 times for each RCP and plot. The initial and climate-induced volume values were classified to volume classes using class limits described above (see also Appendix A). Considering each of the simulation periods separately, it was computed how many times of the 1000 draws the plots with the enhanced growth increased to a higher volume class than with the initial transition probabilities. Multiplying this proportion with the area represented by a plot gave the total area estimate that was expected to transit faster under the given RCP. These area proportions were translated to transition probabilities similar to those based on the initial pairwise NFI observations (Section 2.3.1). In the EFDM simulations, the initial J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 215 transition probabilities (the single probability matrix derived from the NFI observations) were used for computations representing the current climate, and RCP and simulation period specific transitions obtained using the drifter for those representing climate change. 2.4. Harvest scenarios and management activities We projected the development of the carbon stock assuming three different harvesting levels for roundwood: Business-AsUsual (BAU), National Forest Programme (NFP) and NonDeclining Volume (NDV). The BAU and NFP scenarios were based on fixed harvesting levels of 61.75 mill m 3 /a and 76 mill m 3 /a, respectively, which were obtained from the National Forest Strategy of Finland (Ministry of Agriculture and Forestry, 2015)by adjusting the proposed levels to the 95% amount of growing stock considered here. The NDV scenario had an adaptive level of harvests that was determined as the level, which did not decrease the volume of the growing stock. The allocation of harvests to different types of forests was determined according to activity probabilities, which give the proportions of management activities in the data classified according to the factor combinations. Two alternative approaches were tested as the allocation of the future harvests. First, a businessas-usual allocation (A BAU ) was obtained as the proportion of areas with no management, thinning or final felling realized during the most recent five-year period, i.e. based on observations from the permanent plots between the two subsequent inventories. Second, a schoolbook-allocation (A SB ) was obtained as the proportion of areas, which were marked in the field with a need to be harvested within the next five years strictly according to the instructions of forest management (Yrjölä, 2002). Fig. 1 depicts the differences in the two alternative harvest allocations. To fulfill the harvesting objectives, the activity probabilities based on both alternative harvest allocations were iterated to yield a harvesting drain of roundwood corresponding to BAU, NFP, or NDV scenarios (Vauhkonen and Packalen, 2017). The activity probabilities used in all analyses were computed assuming growth rates of the current climate, i.e., using transition probabilities derived from the pairwise NFI observations. ‘‘Uncertainty” in the impact analyses therefore originates from potential climate or management induced additional growth that occurs on top of the growth that is expected based on the transitions observed at the time when making the harvesting decisions. 2.5. Output coefficients Coefficients for the mean values of area classes were determined to quantify the timber assortment drain and carbon stocked in or extracted from the forest. Similar to Vauhkonen and Packalen (2017), we derived the timber assortment drain by computing the relative proportion of logand pulpwood proportions of the entire volume separately for final fellings and thinnings, using all NFI plots or those plots for which the thinnings were simulated, respectively. The biomass in components (stem, branches, foliage, stump, roots) was calculated for each plot (Repola, 2008, 2009). To obtain the carbon content, the biomasses were multiplied by species-specific expansion factors (around 0.5; see Table 1 of Pukkala, 2014). 3. Results 3.1. Effects of climate-induced additional tree growth to the transition probabilities The transition probabilities that were adapted to the changing climate differed between the RCPs as expected. Fig. 2 depicts these differences as the proportion of area that is expected to transit faster than with the probabilities derived from the pairwise NFI observations. Fig. 2 indicates that the area expected to transit faster increased according to the time steps of the simulation, but this proportion varied according to site characteristics such as fertility. When the source and target classes of the faster transitions are examined at the class-level (Fig. 3), it can be seen that average transitions were usually equal to those based on the initial pairwise observations. However, there were more frequently jumps of more than one class than was observed in the case of the transition probabilities derived from the pairwise observations. The results above are based on using RCP4.5 as an example scenario. With respect to other RCPs, the results did not essentially differ except for the magnitudes of climate-induced changes. The conversion of areas to carbon using the output coefficients yielded a development pattern that can be assumed to realistically mimic the carbon dynamics in boreal forest. Among individual carbon (or biomass) components, stem or foliage grew most or least rapidly, respectively. The differences between site types were in the order of 5.7–8 percentage points, species 8.5–11 pp, and geographic areas 8.4–10.7 pp of the component growth, and climate change amplified these figures. However, the main differences Fig. 1. Proportions of management activities in the business-as-usual (A BAU , left) or schoolbook (A SB , right) harvest allocation. Refer to Appendix A for the definition of volume classes. 216 J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 were related to age-size distributions: Fig. 4 shows that climateinduced predictions had similar trends than those derived from the pairwise NFI observations. Especially the growth of the least mature forests was amplified by the climate-adapted transition probabilities. 3.2. Development of the carbon stock and the related uncertainties due to climate and management Figs. 5–6 present the results of the development of carbon stocks and drains, when the transition probabilities described in Section 3.1 were applied together with the activity probabilities for harvests. The difference between Figs. 5 and 6 is that in the former, the activity probabilities (Section 2.4) were iterated in every simulation step to yield the desired harvest goal, i.e., an equal volume of roundwood was harvested in every simulation period. In Fig. 6, the activity probabilities were fixed to the level that yielded the harvest goal in the beginning of the simulations, i.e., an equal proportion of area was always harvested, whereas the volume harvested after the first period depended on the (in-)growth of forest area to the specific class. In all cases, the management decisions were made according to the transition and activity probabilities derived from the original pairwise NFI observations, which is illustrated using grey bars in Figs. 5–6. The lines in the figures depict alternative courses of development, which were obtained by replacing either transition probabilities with those induced by climate (Section 3.1) or business-as-usual harvest allocation with that based on schoolbook (Section 2.4) or applying both of these changes. 3.2.1. Carbon dynamics under current climate and business-as-usual management Based on the business-as-usual transitions and activities, with activity proportions iterated to yield the desired harvest goal (Fig. 5), harvesting 61.75 mill m 3 of roundwood per year (BAU scenario) increased the carbon stock from 810 to 1065 mill tonnes (31%) by the end of the simulation. Harvesting 76 mill m 3 according to the NFP scenario first slightly increased the stock, but ended up to reduce the stock to 757 mill tonnes (7%) in the end of the simulation. The requirement to not decline the volume resulted in a linearly reducing amount of harvests and, subsequently, carbon extracted. The results differed considerably, if proportions of area yielding the desired harvest level were fixed in the beginning of the simulations (Fig. 6). The differences were also most pronounced with the BAU harvesting scenario, where the level of harvests increased by time. This was because of higher proportion of forest that matured due to low initial harvests and the fixed proportion of this area to be harvested according to the activity probabilities. As a result of increasing the harvests, however, the growing stock decreased. With the NFP scenario, the development was fairly similar as described in the previous paragraph. With the NDV scenario, the growing stock decreased unlike when the activity probabilities were continuously iterated during the simulations, whereas the harvest levels reduced less than when iterated. 3.2.2. Level of uncertainties due to climate and management Both the improved growth and the change of harvest allocation from A BAU to A SB increased the harvesting possibilities in all scenarios described above. Compared to business-as-usual, the climate change alone resulted to 5–33% higher carbon stock in aboveand belowground forest biomass in addition to obtaining 8–20% higher harvest drain. The aforementioned levels depended on the harvesting goals and if the activity probabilities were iterated during the simulations (Figs. 5–6). A comparison of the climateand management induced effects shown in Figs. 5–6 indicates that the impacts of climate change were partly related to those produced by the different harvest allocations. Although the same amount of roundwood was harvested in A BAU and A SB , the amount of total carbon (including all components in addition to tree stems) affected was smaller in A SB ;an observation that is further examined in the next section. The comFig. 2. The proportion of area that was expected to transit faster based on the transition probabilities adapted for expected climate (RCP4.5) compared to those based on pairwise NFI observations. The trend lines represent the area proportions obtained using the median value of the climate-induced temperature distribution, while the vertical lines show the variation between the 10% and 90% values of the distribution. The legend refers to dominant species in different site fertility classes – the fifth class including all poorly productive forest land was omitted from the figure. Fig. 3. Transitions from the initial (Class#, t0) to the subsequent volume class (Class#, t1) in the initial pairwise observations (grey symbols) and those adapted for climate-induced growth (black symbols; transitions adapted for RCP4.5 are shown). The thick horizontal lines depict the median, the bottom and top of the boxes the interquartile range between the 25th and 75th percentiles, the whiskers the lowest data within 1.5 times the interquartile range, and circles the data not included in categories above. J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 217 bined effects of climate and management were more pronounced in the case, where the activity probabilities were fixed in the beginning of simulations, for which reason the text below mainly refers to Fig. 6. In it, a corresponding increment in the growing stock could be observed due to climate change (about 8% using the expected values of RCP4.5) as if the harvest allocation was changed from A BAU to A SB (7–9% depending on the harvest target). In addition, if both the climate change and the change in harvest allocation occurred, it increased the amount of carbon harvested to approximately the same level as with climate change alone by the end of the simulation, but with a completely different development pattern during the simulation period (see Fig. 6). 3.2.3. Uncertainties related to the dynamics of individual carbon components A decreasing amount of the total carbon extracted could be observed based on Fig. 5 (e.g., lower left panel), even if the activity probabilities were iterated to yield the desired harvest goals. This observation and also the considerable differences in total carbon stocks and drains depending on the harvest allocation (Fig. 6) might not seem intuitive, when considered as total carbon, but can be better reasoned when the analyses are broken down to the level of individual carbon components (Fig. 7). Although the harvest decisions are based on the amount of roundwood, the total carbon affected by the harvests includes foliage, branches, stump, Fig. 4. The mean annual increment of the total biomass in the NFI plots as a function of age – refer to Appendix A for the definition of the age classes. The thick line depicts the increment in the initial pairwise observations and the thin lines the increments for the eight simulation periods, when the pairwise observations were adapted to the RCP2.6 (left), RCP4.5 (middle), and RCP8.5 (right) climate scenarios. Fig. 5. The development of total carbon during the eight simulation steps, when the activity probabilities were iterated between the individual steps to obtain harvesting goals BAU (left column), NFP (middle column), and NDV (right column). The grey bars depict the development based on business-as-usual harvest allocation (A BAU ) and transition probabilities derived from the pairwise NFI observations. The black ‘‘error bars” are based on the same transition probabilities, but schoolbook-allocation (A SB ) of the harvests. The lines above the bars depict the development, when the transition probabilities are adapted for the climate: the red, black and green lines refer to RCP2.6, RCP4.5, and RCP8.5 scenarios. The solid or broken lines refer to the use of A BAU or A SB , respectively, and there are three lines per climate scenario and harvest allocation, depicting the 10%, 50% and 90% values of the expected temperature distribution. 218 J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 and roots of the trees harvested. The proportion of the different components affected is related to the different allocation of the harvests (Fig. 1): in A SB , most of the drain is obtained from regeneration harvests, which are essentially final fellings of evenaged forests. The carbon harvested in A SB therefore originates from more matured forests, where a higher proportion of the total carbon is stored in the stem wood. Thus, even if both harvest allocations yielded the same level of roundwood, the level of total carbon extracted varied due to the different proportion of carbon components in the forests subject to the harvests. 4. Discussion 4.1. Modeling climate-induced forest growth variations using areabased Markov chain models To account for the effects of climate change, some earlier simulation studies carried out in Finland (Nuutinen et al., 2006; Kallio et al., 2013) have used the models of Matala et al. (2005, 2006) to predict how increasing annual mean temperature and ambient CO 2 affect the forest growth. When incorporating these effects to the forest development scenarios, the studies mentioned above have considered either immediate or gradual increase in the growth, neither of which is realistic according to the growth patterns observed in the past (Pasanen, 1998; Henttonen et al., 2017). Alternative approaches considering the stochasticity of growth have also been presented (Pukkala and Kellomäki, 2012), but additional considerations were needed to implement these effects with matrix models that assume the stationarity of the transitions. Our approach used the RCPs and subsequent climate projections calibrated for Finland (Ruosteenoja et al., 2016) to derive probability distributions for increase in CO 2 and annual mean temperature. The distributions were sampled to provide CO 2 and temperature values to be used as predictors of the growth effect (Matala et al., 2005, 2006). By means of sampling and adding the obtained growth trend with stochastic variation, we were able to account for the uncertainties related to the RCP predictions. Finally, to be applicable in our simulations, the growth series were not used as such, but as classified to respective volume classes. The classification step fundamentally ‘‘smooths” the growth series, as the additional growth modeled for every plot is not transferred to the projections as such. Instead, the area fraction represented by the plot is divided according to the probability of the plot to transit from the original volume class to a higher class. The probability depends on both the distance of the initial volume value to the class limit and the magnitude of the relative scenario effect, which further depends on factors such as site type and tree competition (Matala et al., 2005, 2006). In principle, this type of probabilistic approach could moderate the model-based predictions of climate-induced growth rates, which may otherwise seem overly optimistic (cf., Pukkala, 2017b). Pukkala (2017b) came to this conclusion using an alternative approach, which predicts tree survival rates in addition to straightforward changes in the climateproduction relationship due to the elevated temperature. The mean annual climate-induced increments of stem volume were in the order of 3.5–12.5%, but as much as 30% during the simulation period (Pukkala, 2017b). In our study, the mean annual increment of the total biomass (computed as an average over the age classes in Fig. 4) was 6.8% and was expected to increase to 8.4– 9.4% in the last simulation period depending on which RCP was assumed. Also the results of the development scenarios, when compared in the end of the simulations, are fundamentally in line with Pukkala (2017b). Even though the figures are overall difficult to compare due to different biomass components and computation Fig. 6. The development of total carbon during the eight simulation steps, when the activity probabilities were fixed in the beginning of the simulation. Refer to the caption of Fig. 5 for the interpretation of the image. J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 219 methods (in Pukkala, 2017b, the figures were based on maximizing the net present value with varying interest rates), the approach tested here is considered promising and to warrant tests with other applications than those aiming at Markov chain models. Adapting the transition probabilities to the changed climate slightly changes the concept of the EFDM, which is to a high degree built upon the Markovian property, where the future state depends only on the present state, not on events that preceded it. This property is not affected, but there is a requirement to estimate multiple transition matrices and apply each matrix separately within the given time steps. Conceptually, the changes result to a time-inhomogeneous Markov chain, which is much less applied or even studied than the theory of homogeneous Markov chains. Markov chain models with random or conditional transitions have been tested in other applications than forestry (Shamshad et al., 2005; Meidani and Ghanem, 2013). To date, Liénard and Strigul (2016) are apparently the only ones to describe time-inhomogeneous Markov chains applied to forest projections. Their model was operated at the forest or patch level, which is not directly comparable to our analyses because of larger areas represented by the matrix cells and the use of transitions conditional to management activities in our paper. Therefore, the work on modeling the transition probabilities described here may provide an interesting contribution as an application of timeinhomogeneous Markov chain models with activity-conditional transitions. Fig. 7. The carbon drain presented in the lower row of Fig. 6 broken down to individual (biomass) components. The columns represent the harvesting goals BAU (left), NFP (middle), and NDV (right). The bars contain the proportions of stem, foliage, branches, stump, and roots, respectively, from bottom to top, and separately for final-fellings (dashed bars) and thinnings. The two bars of each simulation period show a comparison of A BAU vs. A SB with transition probabilities derived from the pairwise NFI observations (first row); or initial transition probabilities vs. those assuming the realization of the 90% value of the expected temperature distribution of RCP8.5 under A BAU and A SB (second and third row, respectively). 220 J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 4.2. Limitations of our results and further uncertainties related to projecting future carbon balance The development of the forest resources was projected until around 2050. According to sensitivity analyses based on historical data or comparisons of scenario projections (Nabuurs et al., 2000; Vauhkonen and Packalen, 2017), similar matrix model projections as applied here could be used for periods of up to 50–60 years. Even if the projections could thus have been made for a slightly longer period and also the climate scenarios extend beyond 2050, the trend in many scenarios changes approximately around 2050. Yet, numerical predictions for the GHG concentrations are provided only for years 2050 and 2100 (Meinshausen et al., 2011). Because the projections were composed of sequences of 5-year periods, we feel that excessive assumptions would have been related only to modeling the climate trend beyond 2050, for which reason the projections were not extended further. By ‘uncertainties’ in our analyses, we essentially refer to future development that cannot be predicted using transition probabilities derived from the pairwise observations based on measurements. Our focus was particularly in the effects of climateand management-induced additional net growth to the carbon stocked in aboveand belowground forest biomass. However, both the emphases on net growth and forest biomass also produce limitations towards the interpretability of our results, which is discussed below with respect to earlier literature. 4.2.1. Negative impacts of climate change Climate warming likely affects not only forest growth, but also health: risks of natural disturbances and even calamities increase. If models for occurrence and degree of damages were available, those could be used as additional activity and transition probabilities, respectively, in our model, and therefore also the negative impacts could easily be assessed with respect to the future scenarios. However, even though continent-specific indications on the increase of both biotic and abiotic damages due to climate warming have been presented (e.g., Seidl et al., 2017), there are no numerical estimates available to be used as probabilities. Using a similar Markov chain approach than in this study, the related effects need to modeled indirectly unless the damages are specifically related to the main axes of the matrices (e.g., volume and age). Although this can be done in the EFDM via output coefficients, their use may add further uncertainties to the projections. The use of models developed by Matala et al. (2005, 2006) in the way described in the previous sections allows direct modeling of volume increase as a function of CO 2 , temperature, and forestspecific characteristics. Except for those models, we are not aware of any other climate-adaptive models for Finland that could be integrated to practical simulation systems making use of forest data collected for extensive areas. Liénard and Strigul (2016), who used a Markov chain based approach, found a divergence of ±5% between the development scenarios for mean biomass of hardwood forests of Quebec, Canada, by the beginning of 2090. However, they assumed the increasing CO 2 and temperature to affect more on fire rates than growth enhancements. In Finland, studies integrating future risks in forest projections have considered especially storm-related damages (e.g., Reyer et al., 2017). Also pest attacks or pathogen infections may be expected to increase, but these are more specific in terms of occurrence areas and species, and may therefore require very delicate species-specific modeling (e.g., Nevalainen et al., 2015). Alrahahleh et al. (2016), using a forest ecosystem model that predicts re-generation, growth, and mortality according to temperature sum, tree competition, and soil, nitrogen, and ambient light and CO 2 availability, concluded that climate warming affected the trees in northern and southern Finland indifferently. They elaborated these findings with discussion on joint climate warming, site-specific water holding capacity, and species-specific responses to these phenomena. As mentioned above, our analyses did not account for either biotic or abiotic damages or disturbances except for added growth. We acknowledge that including only positive effects of climate change may be simplistic, but the discussion above also suggests the complexity of considering all possible climate-related impacts. 4.2.2. Other components of carbon balance than biomass The analysis presented in Section 3.2.3 explains the variations in carbon dynamics, when all carbon (or biomass) components are included as the total carbon. Although we acknowledge that the possibility to utilize all these components especially from thinnings can be questioned, we found this analysis beneficial from two aspects. First, an idea on the computational uncertainties involved is provided: the accuracy of biomass models and conversion factors regarding the different components may vary (cf., Neumann et al., 2016). Second, if the components are not extracted and used for biomass-based products, those provide the litter and debris that accumulates as dead organic matter and affects the soil carbon. Our analyses are not complete with respect to the total carbon balance of forests, as we did not explicitly include the carbon sequestered in the forest soil and products. Regarding soil carbon, however, mainly the initialization of the carbon pools is problematic, whereas simulating the decomposition can be based on existing soil carbon models (cf., Pukkala, 2014, 2017a; Akujärvi et al., 2016). However, all aforementioned studies assume that the current climate prevails and no models similar to those applied for living biomass in Section 2.3.2 can be found from the literature for predicting if the decomposition should be assumed to accelerate or slow down in the warming climate. The dynamics of these pools could be estimated by means of coefficients (Pukkala, 2014; Heinonen et al., 2017), but when no climate-adaptive models for the decomposition exist, the changes in the soil carbon would only be related to the varying amount of harvests. However, the effects of harvesting to these stocks can, to a certain degree, be deduced from the earlier studies (Pukkala, 2014, 2017a; ZubizarretaGerendiain et al., 2016). In the simulations of the development of biomass, soil, and product pools under four management scenarios (Pukkala, 2017a), the soil carbon varied much less than carbon stocked in biomass and products. In those simulations, a fast decomposition of all pools was started after a harvest, resulting to a negative total carbon budget in the short term (three to five decades), but a positive budget in the longer term due to the sequestration in the growing trees and harvested products. Thus, according to Pukkala (2017a), the conclusions depend on the time horizon and also on how much weight is set to the substitution effects, i.e., the reduction of consumption of fossil fuels due to wood-based products (see also Pukkala, 2014). The examples above illustrate the complexity and assumptions required to model the carbon balance beyond the living biomass stocks, in changing climate and in the national scale. Due to focusing on the carbon stocked in aboveand belowground biomass, our study cannot be directly compared to those reported earlier. Overall, the multitude of studies and different approaches indicates the challenges in the related modeling task. Earlier carbon balance studies to provide instructions for forest management (Pukkala, 2014, 2017a; Zubizarreta-Gerendiain et al., 2016) were focused on single forest stands or small forest properties (up to around 1000 ha). Heinonen et al. (2017) used NFI11 data and considered the carbon balance of entire Finland, but assumed different harvest allocation and no climate change. The study by Alrahahleh et al. (2016) relied on several assumptions behind the forest ecosystem model used. Compared to that, the J. Vauhkonen, T. Packalen / Ecosystem Services 33 (2018) 213–224 221