Accounting for multiple forcing factors and product substitution 1 enforces the cooling effect of boreal forests 2 3

There is dispute over the climate change mitigating effect of boreal forest management due to the contrasting 20 influence it has on different vectors influencing radiative forcing (RF). For the first time, this study has combined the 21 estimated effects of carbon sequestration in forests and wood products, the surface albedo of forests, the direct and 22 indirect forcing of secondary organic aerosols and the avoidance of fossil emissions by product substitution, both in 23 the current and predicted 2050 climate. The aerosol effect was comparable in magnitude to that of carbon sequestration 24 and increased in importance in a warmer 2050 climate. Harvesting decreased the formation of climate cooling aerosols. 25 The aerosol effect was also larger than the opposing impact of increased surface albedo due to clear cutting in conifer 26 forests. When all above mentioned RF factors were accounted for, the RF of conifer-dominated stands was less 27 negative than that of broadleaf-dominated stands, despite the higher carbon sequestration of the former. Considering 28 also the cooling effect of product substitution, the differences in the RF impact of management alternatives that 29 maintained or increased forest biomass were small. However, the outcome depended heavily on the wood use pattern 30 and the assumed product substitution. A substantial increase in harvest with a clear increase in the share of small 31 dimension fiber and fuel wood use led to a clear climate warming effect in the simulations. 32 33 Biogeosciences Discuss., doi:10.5194/bg-2017-141, 2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c © Author(s) 2017. CC-BY 3.0 License.


Introduction
The boreal biome comprises 1/5 th of the terrestrial carbon sink and about a third of terrestrial carbon stocks (Pan et al. 2011), and its associated effects on surface albedo (A) and secondary organic aerosols (SOA) may have a considerable climate impact (Betts 2000, Tunved et al. 2006).Recently, the cooling impact from afforestation in the boreal zone has been questioned, because models have indicated that the climate warming effect due to increased radiative forcing (RF) from boreal forests' low albedo (Betts 2000) may exceed their cooling impact by carbon sequestration (Avila et al. 2012).However, this argumentation does not account for the possible climate cooling effect of boreal forests via their influence on SOA formation, nor product substitution.Forest ecosystems are producing large amounts of Biogenic Volatile Organic Compounds (BVOC), which leads to SOA production and a greater climatic cooling effect (Tunved et al. 2006) than, e.g., grasslands.Different forest management strategies (e.g.variations in harvest intensity, selection of regenerated species and length of rotation period) can have significant impacts on carbon sink and storage (Pihlainen et al. 2014), as well as on the non-carbon effects, i.e., albedo (Matthies and Valsta 2016) and SOA formation.
Over half of the boreal forest area is currently managed or commercially operable, providing up to 17% of the global industrial roundwood harvest (Burton et al. 2010, Gauthier et al. 2015).Since wood products can be used for substituting more carbon intensive non-wood alternatives, the use of forests provides a potentially very useful and cost-efficient tool for mitigation of greenhouse gas emissions (Sathre and O'Connor 2010).The current harvest level of European forests is substantially lower than their wood increment (Forest Europe…2011), which could be seen as an opportunity to increase biomass utilization to replace fossil fuels.Given the multitude of climate impacts associated with forest management, an assessment of its net influence cannot, therefore, only be reduced to carbon storage or the substitution of fossil fuels by bioenergy.A more comprehensive consideration is required that includes both substitution and non-carbon effects.
In this study, we evaluate four key climate impacts of forest management: (1) carbon sequestration (in trees and soil, i.e. forest ecosystems, and harvested wood products), (2) surface albedo of forest area (A), (3) forest originating aerosols (SOA), and (4) avoided CO2-emissions from wood energy and product substitution (PS).We calculate the net of these effects at both a single stand and regional level, and also consider the effect on global climate if this management is scaled up to the boreal forest level (Fig. 1).Our aim was to have a more comprehensive analysis of climate impacts than the carbon balance.We were particularly interested in on the climatic balance of forest originating SOAs and forest related albedo, and how the avoided CO2 emissions through product substitution compensate the carbon sink reduction caused by wood harvesting.Finland was used as a case study region based on the wide availability of data and validated models for predicting the dynamics of managed boreal forests (Salminen et al. 2005, Mäkelä et al. 2008a,b, Tuomi et al. 2011, The Finnish Statistical…2014).

Overall description of the study
First, we simulated RF for boreal forest management in Finland at the stand level for three dominant boreal tree species in single species stands (Table 1, Scots pine (Pinus sylvestris), Norway spruce (Picea abies), and silver birch (Betula pendula)), and for different site fertility over 100 years under both the current and projected 2050 climate (SRES A2 scenario, ensemble median climate model inmcm3.0,IPCC 2007) (Fig. 1, a and b).Then we estimated RF development at the regional level in Finland, initializing the simulations with recent forest inventory data (The Finnish Statistical… 2014) for both climate conditions (Fig. 1, c).Forest management scenarios were expressed as the ratio of forest harvest to the present current annual stem wood increment (CAI, same absolute harvest levels also in 2050 climate) at the regional level (Fig. 1, d).For stand level simulations, the local impact of a management change was expressed as the change in RF at the top of atmosphere per square meter of forest.The RF values were then compared to a forest just after clear-cutting (bare land) which was used as a reference here (year 0, Fig. 2 and Fig. 4, a and b).
We also calculated cumulative radiative forcing over time (CRF) for illustrating the cumulative effect of separate forcing agents (Fig. 4, c and d).For regional level simulations we reported the change in global RF due to the changes of Finnish forests from their current state (Fig. 5, The Finnish Statistical… 2014).

Estimation of the forest development and the radiative forcing effects of carbon sequestration
Forest growth in Finland under the current climate was simulated with MOTTI, which is an empirical stand-level analysis tool and decision support system for forest management (Salminen et al. 2005).MOTTI was used for calculating tree growth dynamics for pure stands of Norway spruce, Scots pine and silver birch for site fertility types in 3 classes (fertile, medium fertile, infertile).For each combination of tree species and site fertility, the forest management practices spanning a stand's whole rotation in the simulations followed the recommendations for private forest owners in Finland (Table 1, Recommendations for…2006).These recommendations are implemented in MOTTI system and could be selected as an option in simulations.In the simulations, each stand type was artificially regenerated by planting with site specific stand density.Thereafter, the tending of saplings, pre-commercial thinning, commercial thinnings and final harvest was performed according to the recommendations.The timing of the precommercial thinning was based on the site specific tree height while the timing of the commercial thinning depended on the site specific basal area limit defined in the recommendations.Whenever this limit was crossed a thinning from below was performed.In the first commercial thinning, the opening of strip roads was mimicked by removing 18 % of the stand basal area.The treatment intensity was defined on the basis of the basal area in such a way that removal of trees decreased the stand basal area to the recommended level.The timing of final harvest was based on either stand age or stand mean diameter (Fig. 2, Table 1).Final harvest was done in a conventional Finnish way i.e. as a clear cut in which trees are topped, limbed and stems cut to the pre-defined lengths at the harvesting site.Harvesting residues, i.e. tree tops, removed branches and stumps were left at the site.The simulation setup covered about 94% of the upland forest growing sites in Finland (The Finnish Statistical…2014).Litter input into soil was calculated by multiplying simulated biomass compartments and turnover rates in MOTTI with biomass expansion factors (BEF) used in the Finnish greenhouse gas inventory (Official Statistics of Finland 2017).The soil decomposition model Biogeosciences Discuss., doi:10.5194/bg-2017-141, 2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0 License.YASSO07 (Tuomi et al. 2011) was used to simulate soil carbon dynamics with the obtained litter inputs.The carbon stores and decay of harvested wood products were computed from the timber harvest assortment information that MOTTI simulates, and species-specific product distributions and life cycle information.For sawlogs, the allocation of roundwood to industrial products was based on Karjalainen et al. (1994) because that source enabled species-wise computations.Karjalainen et al. (1994) was the most detailed source available and resulted in the breakdown of production given in Table 2. Pulpwood-sized wood was added to fibre products from sawlog residues and consumed for paper and packaging products.Storage of carbon in products was based on four lifecycle categories with exponential decay, given in Karjalainen et al. (1994).
The total uncertainty of modeled forest carbon dynamics was assumed to be 15% in our analysis.This value consists of parametric uncertainty of MOTTI (Salminen et al. 2005) and Yasso (Tuomi et al. 2011) models, and errors in BEFs between different stand age classes.
The marginal annual change of RF due to annual carbon stock changes (in trees, soil and harvested wood products) was estimated following Lohila et al. (2010).The RF is defined here as a local change in the radiative energy balance (W m -2 ) at the tropopause per one square meter in response to the change of the state of 1 m 2 of forest area for the stand-wise simulations.The RF estimates due to the annual changes in biomass were converted to CO2 equivalent emissions/sinks assuming 50% carbon concentration in dry woody biomass.The decay of atmospheric CO2 emissions was estimated using the lifetime function f, (IPPC 2007): where t is time, aj,j = 1, 2, 3 are weights, and τj are time constants (IPCC 2013).From these the cumulative changes in atmospheric CO2 due to different harvest scenarios were estimated and dynamics of RF impact were derived.
To illustrate the cumulative warming impact of forest management within a defined timeframe (T) we calculated the cumulative radiative forcing (CRF) (  ) as: Here   is RF related to effect  (CO2, A, SOA and PS) and  = 100 yr.Accordingly, the total CRF is: For the landscape level simulations we estimated the impact of the management of the whole forestry area of Finland for the whole globe (Sect.2.7).

Estimation of the climate change effects on the forest carbon sequestration
To analyse forest growth and carbon sequestration under climate change, the process-based models OptiPipe (Mäkelä et al. 2008b) and PRELES (Mäkelä et al. 2008a) were used.In the OptiPipe model, tree volume growth was obtained through the allocation of gross primary production (GPP) depending on the C:N ratio, temperature sum and soil nitrogen availability.Soil nitrogen availability increases concurrently with the increasing temperature sum.The potential increase of soil nitrogen availability was set higher in simulations for fertile forest sites than for poorer sites (Mäkelä et al. 2008b).The GPP prediction in PRELES is based on the concept of light use efficiency.The model uses daily values for the absorbed photosynthetically active radiation (aPAR) and modifiers for atmospheric CO2 concentration, light, temperature, vapor pressure, soil water and phenological state for GPP calculation (Peltoniemi et al. 2015).The obtained relative changes in volume growth were used to modify the growth functions of MOTTI for the 2050 climate.
The GPP potential of Finnish forests, used as an input in OptiPipe, was estimated using PRELES for the mean climate In the delta change method, the projected climate changes are added/multiplied (depends on the variable) to the observed climate variables in the reference period.It is also assumed that the climate model bias remains the same in the simulations of future climate (Meehl et al. 2007).Firstly, the monthly long-term changes in air temperature (T, °C), precipitation (R, %), vapor pressure (hpa), global radiation (%) between the reference period  and future period for each GCM and emission scenario were calculated.The GCM grid cell center points were re-projected to the projection of the grid of the observed database and monthly changes were bi-linearly interpolated to estimate values for the center points of the 10 x 10 km grid cells.For each grid cell, monthly changes were linearly interpolated to daily changes, which were added to the observed time-series.

Estimation of the radiative forcing effects of product substitution
The avoided CO2 emissions related to wood utilization were computed separately for harvested sawlogs and pulpwood.For sawnwood, the substitution factors of individual studies listed in Sathre and O'Connor (2010) were used as data points to compute the average values for avoided carbon emissions per carbon in raw material.They were 0.913, 0.905, and 0.819, for Scots pine, Norway spruce and silver birch, respectively.For determining the uncertainty range of the values, standard errors of 0.566, 0.560, 0.507 were used for each species, respectively, calculated from the individual studies reported in Sathre and O'Connor (2010).used as in case of forest carbon accounting (Sect.2.2) for the decay of avoided emissions effect.Also RF due to the avoided emissions was estimated following Lohila et al. (2010).

Estimation of the radiative forcing effects of surface albedo
Models of forest albedo were estimated for an area located in central Finland based on MODIS MCD43A3 blue-sky albedos (Schaaf et al. 2002) and forest resource data produced by the Natural Resources Institute Finland (Tomppo et al. 2008).Regression models were used to estimate the tree species specific forest albedos for different volume thresholds utilizing information on the fractional covers of different forest types within the MODIS pixels (Kuusinen 2014).The land cover data were divided into five components: clear cut (growing stock ≤ 5 m 3 ha -1 ), young stand (pine or spruce forest with growing stock > 5 m 3 ha -1 but < 60 m 3 ha -1 ), pine forest (growing stock ≥ 60 m 3 ha -1 ), spruce forest (growing stock ≥ 60 m 3 ha -1 ) and deciduous broadleaved forest (growing stock > 5 m 3 ha -1 ).Species-specific albedo for Scots pine and Norway spruce were assumed to follow a stepwise function during the total rotation.Albedo was assumed to be highest in open clear cuts and linearly decrease for Scots pine and Norway spruce stands until the growing stock reached 60 m 3 ha -1 .Deciduous broadleaved species (mainly birches) albedo was noted to be insensitive to changes in growing stock, so it was estimated as one value for the total rotation (stand volume > 5 m 3 ha -1 ).Monthly means of component albedos were calculated for February-September, but due to low solar zenith angles during winter, there were no good quality albedo retrievals available from October to January.For December and January, the same albedo values were used as for February; October was given albedo equivalent to that of September and for November albedo was linearly interpolated between the values of October (September) and December (February).
The resulting albedo values were translated into net shortwave radiation at the top of atmosphere using ECHAM5 radiative transfer model.The method is explained in more detail in Matthies et al. (2016).The uncertainty of the albedo impact on RF was estimated to be 25%.This albedo uncertainty includes the differences between vegetation types and in radiative transfer.The uncertainty of the RF from albedo of different forest types was estimated as a proportion of their RF in a clear-cut area (used as a reference here).This assumes that changes in RF are directly proportional to changes in albedo.The analysis is based on published values of RF and their errors (Table 1 in Kuusinen et al. 2013).

Estimation of the radiative forcing effects of forest-originating aerosols
Measurements in a boreal Scots pine forest (SMEAR II, Hakola et al. 2012, Bäck et al. 2012, Aalto et al. 2014) and literature values of BVOC emission potentials for different species (Hakola et al. 1998, 2001, 2003, 2006, Smolander et al. 2014) were used to simulate SOA formation for forests of different age and species composition using the one-  2011), and updates and validations are presented in Mogensen et al. (2015).The emissions of organic vapors from the canopy were calculated by a modified version of MEGAN 2.04 (Guenther et al. 2006, Smolander et al. 2014).The estimated emissions are highly dependent on meteorological factors (in particular temperature and light), forest leaf area and leaf biomass, and furthermore, the emission potentials of individual organic compounds that are specific for individual tree species (Mogensen et al. 2015).The data on Silver birch and Norway spruce emission potentials was much less than those of Scots pine, and therefore simulations include more uncertainty regarding these species.
The SOA effects were calculated for three different ages of pine, spruce and birch forest (Fig. 3 and Table 3) and the values over the stand development were interpolated from those stages.The uncertainty range of SOA RF was adopted from the study by Spracklen et al. (2008).They found a RF range from -1.6 to -6.7 Wm -2 for the SOA effect for the difference between no forest and closed canopy.Our simulations gave an average value for the same difference of -3.8 Wm -2 when the current species distribution in Finland was assumed (50% Pine, 35% Spruce and 15% Birch).This yielded an average range of ± 70% for the SOA effect.The aerosol dynamics module in the SOSAA model was based on the University of Helsinki Multicomponent Aerosol model (UHMA, Korhonen et al. 2004) and described in a recent publication by Zhou et al. (2014).
For the climate change effect in the SOSAA, the daily weather for year 2050 was predicted using the SRES A2 emission scenario and the csiro_mk3_5 climate model (0.34°C lower annual mean temperature in 2050 than projected by ensemble median model).The simulated meteorological parameters in SOSAA were then nudged towards these predictions.The BVOC emission changes were calculated using the temperature dependency as expressed in MEGAN 2.04.Outputs from various climate models were used, in order to constrain SOSAA by the expected ambient concentrations of O3, SO2, CO, NO, NO2 and CH4 in year 2050.The predicted mole fractions of SO2 and O3 were obtained from GFDL-CM3 from the CMIP5 data archive (http://cmip-pcmdi.llnl.gov/cmip5/index.html).The predicted mole fraction of CH4 was taken from CESM1-WACCM (from NCAR).There exists no predicted mole fractions of NO, NO2 and CO for year 2050, therefore the annual mean emission rates of NO and CO from IPCC RCP4.5 were utilized.It was assumed that there is an identical change in the NO2 concentration as in the NO concentration.This means that the 2010 measured concentrations were multiplied with the following factors: The changes in direct RF from aerosols were calculated between each canopy condition of the forest and species (Table 3).The method is based on the supplementary material The monthly mean cloud cover fraction was calculated from ECMWF reanalyzed low cloud cover fraction, which is also expected to have contributed towards the uncertainties.The changes in indirect RF from aerosols were calculated similarly as above between each canopy condition of the forest and tree species (Table 3) based on the method presented in Kurtén et al. (2003).The results in both direct and indirect forcing calculations were averaged for annual values in W m -², which represents the total amount of the radiative effect for 1 square meter of forest.

Estimation of the radiative forcing effects of forest area simulations
The forest area simulations were calculated for the three main boreal species (Scots pine, Norway spruce and silver RF changes of forest area over the simulated time period, i.e. surface area of each species and site type combination was multiplied with stand specific RF change curve (Fig. 4).The forest area RFs were areal integrals of stand RF caused by carbon sequestration, surface albedo, aerosol and product substitution effects expressed relative to the initial state of the forest.No spatial interaction between forest stands was assumed.Under the 2050 climate (Sect.2.3), the forests were initialized using the same site and age distributions as in the current climate.However, the initial volume and growth of these classes corresponded to the new equilibrium in 2050 climate by modifying the growth functions in MOTTI.

Results
Forest management has a clear effect on radiative forcing (RF) regionally.The simulations show that both warming and cooling can be attributed to management options, including the choice of tree species and the rotation period.The RF was more negative with decreasing harvest levels, which is mainly attributable to higher carbon sequestration and SOA formation.Substituting fossil fuels and fossil based material generated accumulation of avoided emissions which reduced greatly differences between harvest schemes.At the stand level, the combination of all effects level resulted in a large negative RF over the rotation (Fig. 4, Table 4).The 100-year average RF was increasingly negative with increasing site fertility, and the contribution of SOA and PS together were over 70% of the total cooling impact in all cases (Table 4).The inclusion of the SOA effect to RF assessment of the biosphere increased the negative RF of boreal forests, i.e. climate cooling influence, and enforced the warming impact of forest harvesting.The cooling impact of PS increased along the larger harvests but depended greatly on the substitution factors (Fig. 4).
In the SOA effect, the increased cloud albedo due to forest VOCs dominated since direct RF effect of aerosol particles ranged only from 3% to 6% in different species.The different effects had different dynamics with respect of forest development.Harvesting caused an increase in RF due to the loss of carbon sequestration and a decrease in SOA, but a decrease due to an increase in PS and also in A in the conifers.The deciduous silver birch stand's A had an opposite effect to that of the conifer stand's, being slightly higher than the open area A, i.e. clear cut had a warming effect in terms of A in silver birch.These differences between spruce and birch stands are illustrated more clearly by the cumulative RF (Fig. 4 c and d).The largest impact in CRF in all stands was caused by PS, due to its cumulative nature.
Under the 2050 climate, the SOA effect over stand rotation was larger than in the current climate (from -2.3 to -5.4 Wm -2 for the herb rich spruce site type) due to enhanced BVOC emissions and subsequent SOA formation, related to increased temperatures.This effect more than compensated for the increase in RF from reduced carbon sequestration (from -3.8 to -2.4 Wm -2 ) due to more rapid biomass turnover (litter and soil carbon, also shorter stand rotation time) in the warmer climate.
Forest management schemes at the regional level were simulated with different annual harvest intensities (50%, 65%, 100% and 130% harvest relative to the present CAI) using the modeled stand level increment curves.Totaling the considered direct effects of forests (immediate A, SOA and CO2, excluding PS), meant that any harvest level below 100% of the current growth of Finnish forests, i.e., harvest schemes that increase the forest biomass from current, had a cooling influence relative to the present state under the current climate (Fig. 5a).The difference between the lowest and highest harvest levels of Finnish forests (50 vs. 130% of CAI) led to a net global RF difference of approximately 0.003 W over 50 years (the time period considered most important for climate mitigation actions) with the two lightest harvest scenarios clearly cooling the climate and the most extreme harvest scenario clearly warming the climate from the present (Fig. 5a).
When the product substitution (PS) effect was also considered, the differences in the RF impact between harvest levels was reduced (Fig. 5b).Particularly, the RFs of the harvest intensities that maintained or increased the forest wood stock were all within the assessed uncertainty ranges.The average substitution factors for saw logs and pulpwood (0.9 and 0.7 kg avoided C emissions per kg of C in wood, respectively) led to comparable RF influences over the first 50 years between PS and the other effects.The higher PS associated with a higher harvest rate largely compensated for the lower direct cooling effect of SOA and CO2 so that the net difference between these harvest levels (50% -100%) was small (Fig. 5b).
The differences in RF between the harvest levels increased further under the 2050 climate.Largely due to changes in difference between the lowest harvest rates was minor (Fig. 6).In the 2050 climate analyses, the PS was not considered since it is a property of the technosphere (comprised of all of the structures that humans have constructed) and thus larger and more rapid changes in it could occur than in the forest related properties of biosphere.Overall, the estimated uncertainties for the non-carbon effects and PS exceeded those for forest carbon balances.

Discussion
Our results show that the rarely considered cooling via the formation of SOA as well as avoided emissions by product substitution (PS) play dominant roles in the climate forcing of boreal forests' management.Together, the climate cooling impact of these two effects was about 70% of the average RF effect over stand rotation, but their response to forest management was mutually opposite; harvesting increased cooling by the PS effect but decreased it by the SOA effect.The analysis showed that a changing climate will lead to an increased importance of SOA.Aerosol and albedo effects were large enough to modify the carbon sequestration based order of dominant species' climate change mitigation potential, shifting the rankings of spruce and birch dominated stands.The CRF curves show explicitly how fundamental the inclusion of SOA and PS effects is for the climatic impact of forest management decisions.The large differences in regional RF values between different scenarios underline the importance of forest management as a driver of climatic forcing.Although the maximum cooling by SOA over rotation was lower than that of cumulative carbon sequestration, the average cooling was comparable since the SOA effect tracked the size of the canopy whereas the carbon effect followed biomass accumulation in trees, litter and harvested wood products.For the coniferous species, the SOA induced cooling exceeded the warming induced by low albedo of forest cover under the present climate.The SOA effect thus increased the climate warming impact of forest harvesting even in the boreal region.The SOA-induced cooling of birch exceeded that of conifers, enforcing the Naudts et al. ( 2016) conclusion, though the leaf biomass is lower in silver birch than in coniferous species.These species-specific differences in RF are mainly caused by their different emission rates of BVOCs (Mentel et al. 2009).The differences increased under the 2050 climate (the net difference of local RF between birch and spruce in herb-rich sites increased from -2.2 to -10.6 Wm -2 without the PS effect) as the relative importance of carbon sequestration decreased due to more rapid carbon turnover, whereas the BVOC emissions and SOA formation simultaneously gained importance.This suggests that tree species selection in favor of broadleaf species with higher BVOC emission rates can be of equal or greater importance than species selection for albedo proposed earlier (Betts 2000, Spracklen et al. 2008, Bright et al. 2014).This result bears significance for the vegetation climate feedback since growth predictions suggest that warming climate favors broadleaf species over conifers (Kellomäki et al. 2008).We would like to point out that the uncertainty related to the Unpublished preliminary results (completed after submission of this manuscript), however, suggest that the monoterpene emission of silver birch used for this study could be in the order of 6-18 times too high.This would decrease the aerosol RF of birch so significantly that conifers forests would be cooler than silver birch forest with respect to aerosol RF.
These results strongly suggest that the SOA effect needs to be included in a full assessment of the climate impacts of forest management, however, the conclusions are conditional to our current understanding of SOA formation.
Numerous atmospheric processes and feedback mechanisms influence the RF of aerosols.For example, it was recently discovered that the composition of aerosol size distribution is influenced by Highly Oxidized organic Molecules (HOMs) (Ehn et al. 2014, Jokinen et al. 2015).Further, the exact effect of clouds on climate depends on cloud altitude, which was not accounted for in our simulations.The average RF values due to SOA impacts are in the middle range of previously estimated values (Spracklen et al. 2008), which is reassuring and gives confidence to the estimates.
Our results suggest a weaker albedo impact than indicated previously (Betts 2000).However, the estimated albedo effect was largest in relative terms in low fertility forests (Forest Resources…2000) which probably correspond better to the Siberian forests where a large albedo impact has been estimated (Betts 2000).Nevertheless, our estimates of the RF due to albedo changes were relatively small.Estimates of the albedo effects of deforestation of boreal forests vary widely in the literature (Betts 2000, Bright et al. 2014, Alkama and Cescatti 2016).Similar to our study, O'Halloran et al. (2012) also derived albedo values from MODIS data for stands of different ages but obtained a somewhat larger albedo difference than our study between forest and no-forest patches.This was mainly due to lower winter and spring albedo of mature forests than in the study applied here (Kuusinen et al. 2012).Our estimates are averages of clear-cut areas and forest stands at different ages for large forest areas and should therefore be representative of prevailing conditions in Finland.Nonetheless, the different studies may reflect differences in forest structure and conditions and empirical data have uncertainties of its' own depending on the length of time series or accuracy of inventory data etc.Therefore, caution must be used when interpreting the results for the whole boreal forest.Overall, the importance of the springtime albedo effect will decrease in a warmer climate with a shorter snow cover period.
The avoided emissions are pivotal for assessing the climate change mitigation efficiency of forest management.
Without it, the lower the harvest the larger the cooling climate impact.Conversely, numerous studies agree that in the long run the avoided emissions of product substitution inevitably lead to climate cooling as they accumulate over time (Sathre and O'Connor 2010).At the time of harvest, carbon stored in the forest is lost and, assuming a given demand of products, fossil fuel emissions are reduced as wood products can act as substitutes.Our results show that in a region of established forest management with sustained harvest without a reduction in the regional wood inventory (Gauthier et al. 2015), a moderate decrease or increase of the harvest rate may not bring about large changes in net climate benefits if the substitution effect is considered in the time range of 50 years.This result is highly sensitive to the Biogeosciences Discuss., doi:10.5194/bg-2017-141,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0 License.substitution factors.Within the time span of our analysis, a tightening climate policy should lead to a lower C intensity of the energy system, leading to lower substitution benefits of biomass (Pingoud et al. 2012).A coefficient of variation ±62% was used here for sawn timber based on the reported values of Sathre and O'Connor's (2010).With high values of substitution factors (i.e., replacement of high emission products with wood) the equal cooling was obtained with low harvesting scenario (50% of CAI) and when harvest equaled growth (100% of CIA), whereas low substitution values required low harvest scenarios for maximum cooling.Extremely intensified forest harvests (130% of CAI), with an associated decrease in harvested wood dimensions led to clear reductions in the associated cooling effect (Fig. 5).This is due to smaller substitution factors of pulp and energy wood relative to those of long-lived timber products.
Hence, the combination of larger harvests and a shift in wood use to products with low substitution factors, such as pulp or bioenergy, is not beneficial from a climate change mitigation viewpoint.
The approximate difference between the climatically 'best' and 'worst' simulated management schemes for Finnish forests without the substitution effect up to year 2050 reached 0.13 % of the global anthropogenic RF in 2011 (0.003 out of 2.29 Wm -2 , IPCC 2013).While this number is low, it is important to note that Finland covers approximately 2% of the total boreal forest area (Burton et al. 2010).Although most boreal forests are currently not managed at the intensity found in the Nordic region, more than half of the biome is commercially utilized (Burton et al. 2010, Gauthier et al. 2015).Accordingly, the influence of the selected management could be around 3% (about -0.08 Wm -2 ) of the current anthropogenic RF in 50 years in the current climate, and even higher with climatic warming.The resulting warming effect with increasing harvests can be viewed as a climate debt that results from forest utilization.It may be counterintuitive that such a debt results even though the harvested forests regenerate.However, regeneration is time dependent and therefore, at the landscape level, forest carbon density and SOA formation decrease due to harvesting, as well as with their associated RF effects.To estimate the net effect, we must also consider how the forest products as a whole influence the climate, in relation to alternative products and production chains.These substitution effects have to be always analyzed against the reference situation.At the stand level analysis, we used clear cut area as the reference situation while mature unharvested forest could be seen preferential especially from viewpoint of the carbon dynamics.This does not, however, affect how different stands compare to each other.The main reason for this selection of the reference was that albedo and VOC emissions of clear cuts are equal between species and rotation times while they differ for mature forests.Also, the modeling of the development of unharvested mature forest is far from not being ambiguous.At the landscape level, we used the scenario corresponding the current harvest level as a baseline (65% of the present CAI), i.e. the continuation of the business-as-usual.When the avoided emissions from product substitution were considered, the differences between the moderate harvest rates (50-100% of current growth) were within the uncertainty range.This would mean that within this range, the choice between more or less intensive forest management neither cools nor warms the climate within a 50 year period, assuming the wood use portfolio does not change.However, the result depends on the prevailing production technologies.While the direct impacts of forest on atmosphere depend on biosphere processes, the substitution effect can change quite drastically with technological innovations, and therefore its long term prediction is difficult.The selection of a reference situation is not a trivial.The results of this study demonstrate that deciding on a climate beneficial forest management requires a holistic consideration of all associated effects, and that appropriate management of boreal forests can have an important climate change mitigation role.The inclusion of the rarely considered SOA effects enforces the view that the lower the harvest the more climatic cooling boreal forests provide.Including the effect of product substitution highlights the importance of the kinds of wood products that are manufactured from harvested wood.Our results, in support of previous similar findings, should act to caution policy makers who are emphasizing the increased utilization of forest biomass for bioenergy over long lasting wood products to combat climate change (Kallio et al. 2013, Söderberg andEckerberg 2013).It has been shown that forest management aimed at maintaining the standing biomass with sustainable timber supply is also preferred by many private forest owners in Finland (Valkeapää and Karppinen 2013).
This study's results demonstrate that such management is also superior to intensive forest biomass use for short-lived wood products in providing the most climate cooling influence.Additionally, sustainable harvest options provide greater consideration for other ecosystem services and an insurance for possible adverse effects of climate change on boreal forest health (Gauthier et al. 2015).Overall, managed boreal forests can provide non-trivial sustained climate cooling without any land-use changes while simultaneously maintaining their other important economical and ecological roles.
model projections of SRES A2 climate scenario(IPCC 2007) (between RCP6.0 and RCP8.5 scenarios,Rogelj et al. 2012) for the year 2050.The constructed climate scenarios combined a daily observed data set on a 10 x 10 km grid over Finland(Venäläinen et al. 2005) with changes in the long-term average, simulated by an ensemble of 8 Global Circulation Models (GCM, CMIP3, Meehl et al. 2007) for the A2 emission scenario for 2011-2040, 2041-2070 and 2071-2100.The development of atmospheric CO2 concentration during 21st century was obtained from the BERN carbon cycle model (IPCC 2013).Details of the construction of climate scenarios are described by Rötter et al. (2013).
dimensional chemical-transport model SOSAA(Boy et al. 2011).The model was constructed to reproduce boundary layer transport, emissions, chemistry and aerosol dynamics.Since SOSAA is a column model, it assumes horizontal homogeneity, with the consequence that simulations were carried out for individual tree stands under the assumption that the gas-and particle phase compounds from one stand did not interact with other stands.The simulated aerosol loading for the atmospheric boundary layer column was then employed to estimate the direct and indirect aerosol induced RF.The indirect effect, which dominates the aerosol radiative effects, is based on the method byKurtén et     Biogeosciences Discuss., doi:10.5194/bg-2017-141,2017   Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0License.al. (2003).The meteorological transport was based on the coupled plant-atmosphere boundary layer model SCADIS(Sogachev et al. 2002, Sogachev and Panferoy 2006, Sogachev et al. 2012).The most important supporting equations are provided in Boy et al. ( of Paasonen et al. (2013) and on the article Lihavainen et al. (2009).Biogeosciences Discuss., doi:10.5194/bg-2017-141,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0 License.
Naudts et al. (2016) concluded that forest management favoring conifers has contributed to climate warming since 1750's, particularly in the European temperate forest.Their conclusion was based on carbon sequestration together with albedo and evapotranspiration impacts.We show here that including the effects of forest management on SOA formation acts in the same direction asNaudts et al. (2016) conclusion.The SOA-induced cooling of closed canopy forest relative to the reference clearcut was comparable to that caused by carbon sequestration over stand rotation.
Biogeosciences Discuss., doi:10.5194/bg-2017-141,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0 License.emission potentials of Norway spruce and especially silver birch is high due to very limited published data.

Figure 1 :
Figure 1: Schematic presentation of the study.We (a) simulated the stand development (changes of carbon storage in trees, soil, and harvested products) of most common Finnish tree species at representative forest types (brown-sub-xeric Scots pine, green-mesic Norway spruce, blue-herb-rich Norway spruce, purple-herb rich Silver birch) at the current and the projected 2050 climate and (b) estimated the local radiative forcing (RF) considering surface albedo, CO2 and secondary organic aerosol (SOA) and avoided emissions from product substitution (PS) effects (yellow arrows indicate direct radiation from sun, red arrows infrared radiation).Using this information and the state of the current forests of Finland (c) we calculated how a change in their state influenced global RF relative to the present assuming different policy options for harvesting (d) (annual harvest 50-130% of current annual forest growth) up to 50 years from the initial state both in the current and 2050 climate.Finally, we estimated our results at the level of managed boreal forests.

Figure 2 :
Figure 2: The dynamics of carbon in different compartments in a) fertile Norway spruce and b) fertile silver birch stands.The clear cut has been done one year before the initialization.Stand development and litter inputs were simulated with MOTTI.Soil carbon was simulated with Yasso07 and was ran to steady-state with the same forest management practices as used in the simulations (Table 1, Recommendations for…2006).From the harvested timber and pulp wood, the portfolio of wood products was formed (Table2) and decay curves calculated on basis ofKarjalainen et al. (1994).

Figure 3 : 816 Figure 4 :
Figure 3: The dynamics of boundary layer aerosols over 100 years in different stands.The dots are the stand ages where the biogenic volatile organic compounds and their effect on the aerosols was modeled (see Sect. 2.6 and Table 3).Lines are linear interpolations between dots.The starting point is bare land one year after clear cut (i.e.final harvest) and the vertical lines depict the year of next clear cut and the second rotation in different stands.

Figure. 5 :
Figure.5: Simulated development of the net global radiative forcing (RF) due to different harvest intensities of Finnish forests (2% of the boreal forest area) relative to their present state in the current climate a) excluding avoided emissions (PS) from product substitution and b) the total impact including PS relative to a reference harvest level at 65% of the present current annual increment (CAI, present harvest level).Different colors represent varying harvest levels given as a percentage of present CAI in Finland in 2013 and the shaded area represents the estimated uncertainty (see text).

Figure. 6 :
Figure.6: Contribution of different factors to the average change of net global radiative forcing (RF) in 50 years due to changes in Finnish forests (2% of the boreal forest area) relative to the harvest level (% current annual increment, CAI).At present, the harvest level is at 65% of CAI.Analyses were done in two different equilibrium climates a) the current and b) the projected 2050 climate.CO2 includes carbon sequestered in trees, soil and forest products, A equals surface albedo, SOA aerosols and avoided emissions from the use of wood products (PS).The PS was not considered in 2050 climate (b) since it is a property of the technosphere and thus larger and more rapid changes in it could occur than in the forest related properties of biosphere.The range bars in a) indicate the uncertainty estimates (see text).

Table 1 : Silvicultural recommendations (Recommendations for…2006) implemented in MOTTI simulation model and 663 applied in business as usual (BAU) scenario simulations. The newly established stands after regeneration on the Herb-rich 664 sites were evenly divided between silver birch and Norway spruce in the simulations while on the other site types it was 665 assumed that stands regenerate with the same species as in previous generation.
Biogeosciences Discuss., doi:10.5194/bg-2017-141,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Table 3 : The characteristics of stands used in modeling the emissions of organic vapors (BVOC) by MEGAN2.04 and the 691 formation of secondary organic aerosols (SOA) by SOSAA. We had standard emission potentials of BVOCs that are 692 different from zero only for April through September for birch due to the leaves fall off the trees and for May through 693 September for Norway spruce.
Biogeosciences Discuss., doi:10.5194/bg-2017-141,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Table 4 : The 100-year average local radiative forcing (RF, W m -2 ) in the current climate of the three most common forest site types in Finland by separate influences (CO2 sequestration in forest stands and wood products (CO2), surface albedo (A), aerosols (SOA) and avoided emissions from the use of wood products (PS) and in decreasing order of site productivity. Values are relative to bare land.
Biogeosciences Discuss., doi:10.5194/bg-2017-141,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 May 2017 c Author(s) 2017.CC-BY 3.0 License.