Bridging mapping and simulation modelling in the ecosystem service assessments of boreal forests:effects of bioenergy production on carbon dynamics

Anu Akuj?rvi ,Anna Repo,Altti M.Akuj?rvi and Jari Liski

Abstract

Keywords: Carbon budget, Ecosystem services, Forest management, Forest bioenergy,Landscape level, Mapping,Modelling

Background

Bioenergy, produced from harvest residues such as branches, tree tops and stumps, is an increasing form of utilising boreal and northern temperate forests (Diaz-Yanez et al. 2013; Scarlat et al. 2015). Bioenergy reduces fossil carbon emissions to the atmosphere. However, increased biomass harvesting reduces the carbon stocks of forest which may partly reduce the climate benefits of forest bioenergy (Schlamadinger et al. 1995; Repo et al.2011; Schulze et al. 2012; Zanchi et al. 2012; M?kip?? et al. 2014). In addition, the increased biomass extraction in the boreal regions has also raised concerns about the degradation and loss of habitats of deadwood dependent species (Bouget et al. 2012). The volume of deadwood is strongly correlated with the richness of threatened species: birds (Virkkala 2016), insects (Martikainen et al. 1999), and fungi(Penttil? et al.2006), which makes it a good biodiversity indicator.

The climate impacts of forest bioenergy have been studied extensively in recent years using simulation models (Wihersaari 2005; Eriksson et al. 2007; Melin et al. 2010; Kilpel?inen et al. 2011; Repo et al. 2011,2012). One approach has been to scale up stand level estimates of the CO2emissions resulting from carbon stock changes to the national level by assuming a uniform age class distribution of the forest stands (Cherubini et al. 2013). This approach approximates the total CO2emissions from a regularly managed forested landscape. The age class distribution has, however, a significant effect on the net emission estimates from energy wood use (Routa et al. 2012). Ignoring the high variability of forest structure and the irregular occurrence of harvests might add inaccuracy to the landscape level estimates of carbon budget. The spatial resolution of the carbon stock changes could be improved by coupling remote sensing-and inventory-based observations of forest characteristics with simulation modelling (Paulick et al.2017). This kind of approach could be applied to illustrate the climate effects of alternative bioenergy production scenarios across an actual landscape where decisions are made. Mapping can reveal the most suitable areas for bioenergy production in terms of resource availability (Verkerk et al. 2019). It could also be applied to identify spatial trade-offs and synergies between bioenergy production and other ecosystem service indicators,such as carbon sequestration and deadwood production(Sacchelli et al. 2013).

The present status and past changes of ecosystem services can be mapped using remote sensing (Vauhkonen 2018; Li et al. 2019) or readily available land cover and land use data (Lautenbach et al. 2011;Mononen et al. 2017). The provisioning potential of ecosystem services is often quantified using simple land cover -based proxies if direct mapping is not possible (Smith et al. 2006; Nelson et al. 2009; Burkhard et al. 2012; Maes et al. 2012). However, the use of land cover -based proxies may simplify the spatiotemporal variability of climate regulation by assuming,for instance, constant carbon stocks per land cover class. This makes the maps prone to errors (Eigenbrod et al. 2010). Empirical and process-based modelling includes more detailed description of forest growth based on measurements or ecological theory.Simulation modelling, in combination with remote sensing data, has been applied to map regulating and provisioning ecosystem services at various spatial scales (e.g. Sitch et al. 2003; Schr?ter et al. 2005;Mina et al. 2017; Holmberg et al. 2019; Verkerk et al.2019). However, the complex structure and data requirements of models might limit their use in ecosystem service assessments at a fine spatial resolution due to the lack of suitable remote sensing data(Lavorel et al. 2017). This indicates a disparity between the scientific knowledge about carbon dynamics and its implementation to mapping tools. Consequently, simple dynamic tools are needed to evaluate the effects of alternative forest management practices on ecosystem services at varying spatial scales (Crossman et al. 2013).

To contribute to bridging the gap between proxybased mapping and detailed simulation modelling in studies investigating the future provisioning of ecosystem services, a framework for quantifying the carbon budget of boreal forested landscapes was developed(Akuj?rvi et al. 2016). It coupled simulated timeseries of carbon stocks with extensive, publicly available forest inventory data (Tomppo et al. 2014). This relatively simple method enabled reliable mapping of the current status of carbon budget at the landscape level, complementing previous approaches. However,the suitability of this framework for investigating the effects of alternative forest management scenarios on the availability of ecosystem services in the future remains to be tested.

The first objective of this study was to apply the previously developed framework for simulating the future development of ecosystem services related to carbon cycle,particularly in the case of bioenergy production. The second objective was to evaluate the suitability of this approach for assessing ecosystem services at the landscape level. The studied ecosystem service indicators were the carbon stocks and fluxes of forest biomass and soil, timber and energy-wood production and the annual coarse woody litter production,used as a proxy for deadwood abundance.They were simulated for a boreal forested catchment in southern Finland for the period 2012-2100.The validity of the mapping framework was evaluated by comparing the simulated estimates with measurement-based data.

Methods

Study area

The study area Vanajavesi catchment in the southern boreal forest zone belongs to the Kokem?enjoki river basin in southern Finland (Fig. 1). It consists of 10 subcatchments of the second level of the Finnish watershed division. Most of the catchment is within the region of Kanta-H?me. The annual mean temperature was 4.2°C and the annual precipitation 637 mm during 1970-2012.The total area of the catchment is 2700 km2, of which 1425 km2is managed forest covered by the model simulations of this study. The proportion of peatlands in Kanta-H?me is 17%, of which nearly all are forested and about three quarters drained (Natural Resources Institute Finland 2020). Despite the relatively high proportion of peatland forests in the study area they were excluded from the simulations. This was because the litter and soil carbon model Yasso operates only on mineral soils. The Natura 2000 protected habitats (Evans 2012) and nature conservation and wilderness areas,protected by the Nature Conservation Act, cover altogether 3% of the catchment area. The protected forests were considered to affect the total carbon budget of the study area little due to their small cover. They were thus excluded from the analyses. Fertile site types, usually dominated by spruce and deciduous trees, covered as much as 90% of the forest area on the mineral soils in H?me-Uusimaa in 2009-2013. The average cover of these site types is about 75% in Southern Finland. Planting or natural regeneration, regular thinning and clearcutting are common forest management practices in the study area. In 2012, about 16% of the harvest removal consisted of energy-wood, of which 30% was spruce(Natural Resources Institute Finland 2020).

Mapping framework

A previously developed framework (Akuj?rvi et al. 2016)was refined to simulate key ecosystem service indicators related to the carbon cycle of forests: carbon stocks and fluxes of biomass, litter and soil, timber and energywood production, and coarse woody litter production.The development of these variables was simulated over the study area for 2012-2100 with standard forest management practices combined with varying levels of harvest residue extraction for bioenergy production. The variability of actual management interventions was considered by stochastic simulations of the stand age class distribution over the simulation period. This could improve the ecosystem service estimates in comparison with some previous studies applying fixed age class distributions (Routa et al. 2012; Frank et al. 2015; Pang et al. 2017).

Fig.1 Land use in the southern boreal Vanajavesi catchment area in Finland in 2012.The location for the fine-scale visualisation of the results is marked with a purple line.CORINE land cover data was provided by Finnish Environment Institute

The mapping framework coupled extensive multisource forest inventory data with simulation modelling of forest carbon dynamics at the stand level. The carbon stock of biomass, timber and energy-wood production and the production of coarse woody litter were simulated using the MOTTI v. 3.3 stand simulator (Hynynen et al. 2002; Salminen et al. 2005). Coarse woody litter production was used as proxy for deadwood abundance and for this purpose it was defined to consist of tree tops and stumps, not including coarse roots. It was coupled with the litter and soil carbon model Yasso15 to estimate the belowground carbon stocks and the changes in them (Tuomi et al. 2009, 2011a, 2011b; Repo et al.2016). The model simulations were carried out separately for altogether 18 forest site type and main tree species combinations present in the study area (Table 1)and generalised for the entire study area using look-up tables and multi-source forest inventory data.

The performance of the mapping framework has been evaluated previously at a fine spatial resolution in southern Finland by Akuj?rvi et al. (2016). According to their results, the method produced reliable estimates of the current status of forest carbon stocks and changes. The previous validity test is a sound basis for applying the framework for scenario analysis in the present study because the model simulations and initial land cover data were identical. However, the reliability of the simulated ecosystem service indicators was also evaluated broadly in this study. The simulated mean estimates for the study area in the beginning of the simulation period in 2012 were compared with inventory-based estimates taken previously from the same area.

Forest management scenarios

In the simulations, forests were assumed to be managed according to the Finnish good practice guidance for forestry with adaptations related to the rotation length. The stand regeneration method, and the timing and intensity of thinning followed the guidelines (Sved and Koistinen 2015). After planting or natural regeneration, the stands were thinned two to three times depending on the mean rotation length, which varied between 65 and 120 years(Table 1). The mean rotation lengths were determined based on the inventory-based range of mean stand age in the study area. They were similar to the upper rangeof the recommended rotation lengths (Sved and Koistinen 2015). Consequently, final felling was conducted only on mature stands. To account for the deviations of practical forest management from the recommended rotation lengths, stochastic variation of the regeneration age was introduced.

Table 1 The mean rotation length assumed for the site type and tree species-combinations present in the study area.The site types were classified according to Cajander(1949)

The effects of forest bioenergy production on the carbon budget were assessed by comparing the extraction of forest harvest residues to a scenario under which they were left on site to decompose. Harvest residues were extracted only from the fertile Norway spruce sites(OMaT, OMT, MT and VT site classes, see Table 1),mimicking the national good practice guidelines for energy-wood harvesting in Finland (Koistinen et al.2016). Three scenarios were studied. In the reference scenario 1, all forest harvest residues were left on site to decompose and no bioenergy was produced. In scenario 2, 70% of the aboveground harvest residues (branches and tree tops) were extracted for bioenergy production from the final felling sites. In scenario 3, 70% of the aboveground residues, stumps and roots were extracted.

Forest stand simulations

The carbon stock of biomass was estimated over the rotation period using the MOTTI v.3.3 stand simulator(Salminen et al. 2005; Hynynen et al. 2014). It is based on empirical growth and yield models describing the structure, growth and management of the most typical site types and tree species in Finland (Hynynen et al.2002). They have been compiled and validated based on extensive forest inventories and field experiments(Matala et al. 2003). Because the MOTTI v. 3.3 simulator uses a time-step of five years, the intermediate annual values were produced by linear interpolation. The annual estimates of timber, energy-wood and coarse woody litter production were produced by converting the dry biomass estimates to fresh volume using a conversion coefficient of 400 kg?m?3, which is a mean for pine, spruce and birch (Alakangas et al. 2016). Coarse woody litter production was assumed to consist of tree tops and stumps, which represented the large-diameter fractions of deadwood forming litter. In Finland, the diameter of tree tops used for bioenergy production is 4-6 cm (Karttunen et al. 2016). In the output of the MOTTI simulations, the basal-area weighed mean diameter of the mature stands varied between 21 and 47 cm,depending on site type. The diameters of different biomass compartments, such as stumps, were not reported explicitly in the model output. In managed forests, a significant proportion of the deadwood is formed during thinning and final felling from forest harvest residues(Er?j?? et al. 2010). In this study, the annual production of coarse woody litter was assumed to indicate the potential deadwood biomass in the future. Deadwood is a widely used indicator for forest biodiversity (Gao et al.2015).

The output of MOTTI v.3.3 was used as input to Yasso15, an improved version of the Yasso litter and soil carbon model (Tuomi et al. 2009, 2011a, 2011b).Yasso15 was used to estimate the soil organic carbon stock, the annual changes in it and the heterotrophic respiration over the rotation period. The carbon input to soil consisted of natural mortality, forest harvest residues and the annual litter production of living trees and ground vegetation. They were estimated with the same method as in the national greenhouse gas inventory of Finland (Ortiz et al. 2013; Siev?nen et al. 2014; Statistics Finland 2018). The annual litter production of the living trees was estimated by multiplying the biomass compartments of standing trees, derived from the output of MOTTI v. 3.3, with compartment- and species-specific turnover rates (Liski et al. 2006). The litter production of ground vegetation was equal to the estimates used in the national greenhouse gas inventory (Muukkonen and M?kip?? 2006). The carbon content of biomass was assumed to be 50%.

In the Yasso15 model, soil organic carbon is divided into four chemical compound groups: ethanol soluble(denoted with E), water soluble (W), acid hydrolysable(A), non-soluble (N). The decomposition rate of each compound group depends on temperature and precipitation and results in the formation of more recalcitrant humus (H). The decomposition rate of woody litter depends also on its physical diameter (Tuomi et al. 2011a).The chemical quality of non-woody and woody litter was derived from previous studies (Ortiz et al. 2013; Siev?nen et al. 2014). A mean diameter of 2 cm was assumed for branches and roots, and 15 cm for stems and stumps. The litter and soil carbon stock was initialised by first running the model to a steady state with average climate and litter production values over one rotation period. The model was then run for a second rotation before the start of the actual simulation period in 2012.Consequently, the litter and soil carbon stock in 2012 demonstrated the situation right after final felling in the bioenergy scenarios 1-3. The average climate in 1970-2012 (Finnish Meterological Institute 2020) and average litter production over the rotation period were used.The mean annual precipitation, temperature, and temperature amplitude were calculated based on the daily observations from the Lammi weather station located in the study area.

To demonstrate the fine-scale spatial pattern of the carbon fluxes, the net ecosystem production (NEP) of forest was visualised for an 8 km2sized subset of the study area (Fig. 1) for 2012, 2050, and 2080 which demonstrate well the temporal variation in this variable. The NEP represents the net uptake of carbon of forest before subtracting harvest removals. It is calculated as the difference between net primary production (NPP) and heterotrophic respiration (Rh). NPP consists of the annual change in the carbon stock of biomass, litter production,harvest removals and natural mortality. The net carbon balance of the forest, i.e. the net biome production(NBP), is determined by the revenues by NEP and expenditures by harvests (see for example Liski et al.2006).

Forest landscape simulation

The simulated stand-level estimates of ecosystem service indicators were connected to spatially explicit information on forest site type, tree species, biomass and stand age. These data were extracted from the Multisource National Forest Inventory (hereafter MS-NFI) dataset representing year 2011 for the studied catchment. The Finnish MS-NFI produces a wide suit of regularly updated forest variables applicable for ecosystem service assessments from landscape up to national scale (Vauhkonen and Ruotsalainen 2017; Kangas et al. 2019). The MS-NFI forest resource maps are based on extensive NFI field plot measurements, high-resolution satellite images and digital maps, and the non-parametric k Nearest Neighbours estimation (Katila and Tomppo 2001; Tomppo et al. 2008a; Tomppo et al. 2008b). The spatial resolution of the MS-NFI data was 20 m×20 m.

To set the initial state of the simulated forests in 2011,each grid cell of the stand age layer of the MS-NFI data was classified based on the forest site type and main tree species present in that cell. The main tree species was determined as the species having the maximum biomass among Scots pine (Pinus sylvestris L.), Norway spruce(Picea abies (L.) H. Karst) and deciduous species, comprising mainly of Silver birch (Betula pendula Roth) and Downy birch (Betula pubescens Ehrh.). Only one site type and tree species combination was assigned to each grid cell. The stand age distribution at the initial state,based on the multisource forest inventory, was not altered in the classification procedure (for more details,see Akuj?rvi et al. 2016).

To project the development of ecosystem service indicators at the landscape level, the stand age layer was updated annually over the simulation period 2012-2100.The site type and main tree species composition were assumed to remain the same as in the initial state. To mimic the variation of the regeneration age in practical forest management, probabilities of final felling were applied for each site type and tree species combination.This approach was chosen because the guidelines are not usually followed strictly in practical forest management. The probabilities were determined by assuming each combination a normally distributed rotation length with a mean following the guidelines (Sved and Koistinen 2015) and with a standard deviation of 10 years (Table 1). This deviation was considered reasonable in comparison with typical rotation lengths in southern Finland. For example, the mean rotation length of a spruce stand growing on Oxalis-Myrtillus site type was 90 years, ranging from 80 to 100 years. In other words, the simulated rotation length within this class varied stochastically among grid cells. On every simulation year, the forest in each grid cell was either spared or felled based on the probability of final felling that year. This procedure resulted in a new stand age layer representing each simulation year. The simulation method produced regular harvesting and a steady behaviour of the stand age distribution yet accounting for the stochasticity of actual management interventions (Fig.S1). To map the simulated ecosystem service indicators,they were joined to the annual stand age layers using look-up tables.

Model performance

To evaluate the performance of the mapping framework,the simulated estimates in the beginning of the study period in 2012-2016 were compared with measurements derived from literature and forest statistics. The litter and soil carbon stock (Rantakari et al. 2012) was measured from permanent sampling plots in southern Finland for the first time in 1985 and again in 2006(M?kip?? and Heikkinen 2003). The carbon stock of biomass, and its change were calculated based on the mean sample plot estimates of the national forest inventory in the region of Kanta-H?me in 2009-2018. Timber and energy-wood removal were measured from H?me-Uusimaa forestry region in 2009-2014. The simulated results represent a sample of forests located within these regions.

The biomass carbon stock change was estimated based on the measured growth and removal of stems. This approach was chosen because in managed forests, the biomass carbon stock change mainly depends on harvesting, and measurements of stem biomass, growth and removal are readily available from the study area.The growth rate of stem biomass was assumed to be the same as that of stem volume. The biomass carbon stock change was then estimated by subtracting the harvest removal (consisting of timber and energy-wood) from stem growth. The standard deviation of the simulated and measured estimates was calculated between the simulation years and the available measurement periods,respectively. The tree data were derived from the statistical database of the Natural Resources Institute Finland(Natural Resources Institute Finland 2020).

The estimates of annual coarse woody litter production could not be compared with inventory-based estimates of deadwood volume because of methodological differences. The simulated estimate represents the potential input or flow of coarse woody litter to the ecosystem. The inventory-based deadwood volume is a stock.Moreover, only snags and logs over 10 cm thick and 1.3 m long are measured. Deadwood, litter and soil organic carbon (SOC) pools are not separated in the output of the Yasso15 model (Tuomi et al. 2011b; Didion et al.2016). A proper comparison of the simulated and inventory-based estimates of deadwood would require modelling the decay of large-diameter deadwood separately from the decomposition of carbon pools (Herrmann et al. 2015). This could be a subject of further development of this research. To illustrate the role of the initial age class distribution on the regional carbon stock estimates, the simulated and measured age class distribution in the beginning of the simulation period were compared (Fig. S2). The simulated estimate was based on the combination of 1)the multisource NFI data from 2011, which is already generalized data based on the field samples and remote sensing data and 2) the applied forest management scenarios over five first years.The measured estimate is based only on the NFI field plot samples.

Results

Carbon stocks and fluxes

The annual mean estimates of carbon stocks and fluxes were presented as means over the studied landscape to illustrate their temporal variation (Figs. 2, 4). The simulated carbon stock of biomass fluctuated between 5.4 and 7.3 kg?m?2over the simulation period 2012-2100,independent of the bioenergy scenario studied (Fig. 2a).In the beginning of the simulation period in 2012, the simulated carbon stock of biomass was about 25%higher than the measured mean in the surrounding region(Table 2). The carbon stock of biomass decreased at a mean rate of ?0.003 kg?m?2?year?1in 2012-2100, independent of the bioenergy scenario studied. The rate of biomass carbon stock change varied between ?0.07 and 0.07 kg?m?2?year?1(Fig. 2b). In 2012, the simulated biomass carbon stock change was almost twice as much as the measurement-based mean (Table 2).

The litter and soil carbon stock remained relatively stable over the simulation period, varying between 8.5 and 8.8 kg?m?2(Fig. 2c). The simulated estimates were within the upper end of the range of the measured litter and soil carbon stock estimates (Table 2). Forest bioenergy production caused a mild decrease in the litter and soil carbon stock compared with stem-only harvest.In 2100, the litter and soil carbon stock was 1% and 4%lower in scenarios 2 and 3, respectively, than in the reference scenario 1. The litter and soil carbon stock increased at a mean rate of 0.007 kg?m?2?year?1, with a range of ?0.003 ?0.017 kg?m?2?year?1(Fig. 2d). The simulated estimate of the litter and soil carbon stock change was very similar to the measured mean (Table 2).

Fig.2 The mean simulated carbon stocks of biomass,litter and soil and the annual changes in these stocks in the study area over the simulation period 2012-2100.The total carbon stock of forest is the sum of biomass,litter and soil carbon stocks.Positive values of carbon stock changeindicate a sink of carbon,and negative values a source. The bioenergy scenarios 1-3 represent varying level of harvest residue removal: 1) noextraction of harvest residues,2)extraction of branches and tree tops for bioenergy production, 3)extraction of branches,tree tops and stumproot systems for bioenergy production

The more biomass was extracted for bioenergy production, the slower was the accumulation of soil carbon.The soil carbon sink was 18% and 59% lower in scenarios 2 and 3, respectively, than in scenario 1 in the beginning of the simulation period in 2012 (Fig. 2d). The differences between the scenarios levelled off towards the end of the simulation period. This was because the residues from the earlier harvests had mostly decomposed and the volume of energy-wood harvests started to decline (Fig. 5b). The changes in the total carbon stock of forest were mainly driven by the changes in the carbon stock of biomass rather than the bioenergy scenarios. The total carbon sink was 3% and 9% lower in scenarios 2 and 3, respectively, than in scenario 1 in the beginning of the simulation period in 2012 (Fig. 2f).

Maps of net ecosystem production (NEP) demonstrate the fine-scale spatial variation of carbon fluxes across asubset of the studied landscape (Fig. 3). The mean annual NEP decreased at first as a result of regular harvesting and, consequently, the increasing proportion of regeneration stands. The forest acted as a source of carbon to the atmosphere in the middle decades of the simulation period because the harvest removals exceeded NEP (Fig. 4). The spatial patches of NEP became more fragmented as more grid cells were harvested(Fig. 3).

Table 2 The simulated mean ecosystem service indicators for the study area in 2012-2016 compared with measurements taken from the surrounding regions(SF stands for Southern Finland, KH for Kanta-H?me and HU for H?me-Uusimaa, respectively). The litter and soil carbon estimates from the BioSoil(BS)data were derived from Rantakari et al. (2012)and those of the National Forest Inventory(NFI) from the statistical database of the Natural Resources Institute Finland (2020)

The spatiotemporal pattern of carbon uptake reflected the changes in stand age class distribution resulting from regular harvesting which was simulated with a stochastic

Fig.3 A fine-scale visualisation of the simulated net ecosystem production(NEP)of forest for a subset of the study area in scenario 1(no extraction of harvest residues for bioenergy production). The results represent the potential annual carbon uptake of forest before subtracting harvest removals.Positive values indicate that the ecosystem is a sink of carbon and negative values that it is a source

Fig.4 The mean simulated carbon fluxes in the Vanajavesi catchment area over the simulation period 2012-2100 with bioenergy scenarios 1-3.NPP stands for net primary production, Rh for heterotrophic respiration, NEP for net ecosystem production and harvest for the timber and energy-wood removal.The bioenergy scenarios 1-3 represent varying level of harvest residue removal:1)no extraction of harvest residues,2)extraction of branches and tree tops for bioenergy production,3) extraction of branches,tree tops and stump-root systems for bioenergy production

RhRhRhalgorithm. As a result of continued harvesting, the proportion of fast growing, middle-age stands declined leading to reduced NPP and, consequently, reduced accumulation of the biomass and soil carbon stocks. The carbon uptake of forest recovered along with the increasing proportion of fast-growing intermediate-aged stands (Fig. S1). The simulated and measured age class distribution differed to some extent in the beginning of the simulation period because they were based on different methods (Fig. S2). The measurements showed a higher proportion of very young and old stands and a lower proportion of intermediate aged stands than the simulations.

Wood production and coarse woody litter

The provisioning ecosystem services were shown as annual sums for the studied landscape (Fig. 5). The simulated mean annual timber production from the thinning and final felling sites was 0.91 mill. m3?year?1in 2012-2100, with a range of 0.73-1.1 mill. m3?year?1(Fig. 5a).The simulated estimate of timber harvest in 2012 was quite similar to the measured mean (Table 2). Both timber and energy-wood production peaked in the late 2050s, as more stands reached maturity, and decreased thereafter (Fig. 5a, b). In scenarios 2 and 3, the annual energy-wood potential from the final felling sites varied between 0.018 and 0.052 mill. m3and 0.052 and 0.15 mill. m3, respectively (Fig. 5b). In other words, the extraction of stumps multiplied the energy-wood potential nearly three-fold compared with the extraction of only branches and tree tops. The simulated estimates of energy-wood removal in 2012 were comparable with the reported mean (Table 2).

The annual production of coarse woody litter from the thinning and final felling sites remained at a stable level in the study area in 2012-2100 (Fig. 5c). It followed the development of the total harvest removal in the study area because harvest residues were extracted only from the fertile spruce sites (Fig. 5a). In scenarios 1 and 2, the annual mean of coarse woody litter production was 0.18 mill. m3, with a range of 0.16-0.20 mill. m3. These two scenarios produced the same estimates because only tree tops and stumps were classified as coarse woody litter and the simulated amount of tree tops was negligible.The extraction of stumps in scenario 3 reduced the annual production of coarse woody litter on average by 4.6% (Fig. 5c). The simulated estimates correspond to the mean annual input of 0.23-0.25 t C?ha?1?year?1depending on the scenario.

Discussion

Effects of bioenergy production on ecosystem services

Fig.5 The simulated mean annual production(in fresh volume)of timber a,energy-wood b and coarse deadwood c in the study area over the simulation period 2012-2100.The bioenergy scenarios 1-3 represent varying level of harvest residue removal:1)no extraction of harvest residues,2)extraction of branches and tree tops for bioenergy production,3)extraction of branches,tree tops and stump-root systems for bioenergy production

A refined framework to simulate ecosystem services related to carbon cycle was applied across a boreal forested landscape. Regular harvesting was the primary driver of the spatiotemporal variation of the ecosystem service estimates because it affected the forest age class distribution. The forests in the studied landscape acted first as a sink of carbon. However, they started to release carbon to the atmosphere when the harvest removals exceeded growth. The net carbon balance turned negative despite the positive values of NEP. The carbon sink recovered along with the increasing proportion of fast-growing intermediate-aged stands towards the end of the simulation period. A similar effect of forest age on productivity has also been observed in previous measurement (Pregitzer and Euskirchen 2004) and modelling experiments (Kohlmaier et al. 1995; ?upek et al. 2010),supporting the findings of this study. It is noteworthy that at a national scale, Finnish forests act as a sink of carbon (Statistics Finland 2018). However, intensive harvests may turn a region to a temporary source of carbon(Pussinen et al. 2009), like shown in this study.

Based on the simulations, harvest residue extraction caused trade-offs between climate regulation, energywood production and habitat provisioning for deadwood dependent species. Extracting forest harvest residues in the final felling caused carbon loss from litter and soil compared with stem-only harvest. The level of harvest residue extraction increased towards the middle of the simulation period as more stands reached maturity. It declined again to the original level when the stand age class distribution approached the initial status. The annual coarse woody litter input to soil declined as a result of stump extraction in scenario 3, suggesting that continuing energy-wood harvesting reduces coarse deadwood formation in the long term, supporting the finding of Repo et al. (2020). Aboveground residue extraction in scenario 2 had a small impact on the simulated coarse woody litter production because the proportion of tree tops was so small.

Both fine and coarse deadwood are important for wood-inhabiting fungi (Juutilainen et al. 2014) and beetles (Jonsell et al. 2007), stressing the importance of retaining both fractions in the forest. In addition, about a quarter of all forest species in Fennoscandia is dependent on the availability of deadwood (Siitonen 2001). According to field experiments, the removal of harvest residues accelerates the breakage and disappearance of coarse logs (Rabinowitsch-Jokinen and Vanha-Majamaa 2010; Herrmann et al. 2015). Therefore, the actual impacts of harvest residue extraction on both carbon sequestration and biodiversity could be more severe than estimated in this study. The extraction of harvest residues and stumps, together with the mechanical breakage reduced the amount of coarse deadwood on average 1-3.7 m3?ha?1(Repo et al. 2020).

Evaluation of the mapping framework

The validity of the mapping framework was evaluated by comparing the simulated and measurement-based ecosystem service indicators in the beginning of the simulation period. Based on the results, the simulated estimates of the biomass and soil carbon stocks were generally higher than the measurement-based estimates.The simulated biomass carbon stock change was, however, close to a previous model-based estimate (Forsius et al. 2016). The results may reflect a deviation of the studied sample from the larger population of forests which the measurements were taken from. The simulations produced the annual change in the total biomass of trees and ground vegetation whereas the measurement-based estimate included only tree growth.Moreover, the simulated estimate of the litter and soil carbon stock included also coarse woody litter unlike the observations, partly explaining the difference.

The deviations between the simulated and measured estimates of carbon stocks probably reflect differences of the simulated and actual forest management. In the model simulations, natural disturbances were absent,thinning, final felling and energy-wood harvest were assumed to occur always on time and only mature stands were regenerated. As a result of these assumptions, the estimates of stand growth, litter production and biomass potential might have been overly optimistic as also discussed earlier by Akuj?rvi et al. (2016). The simulated soil carbon stock change was very similar to the measured mean. The simulated timber and energy-wood removal were also in the order of magnitude with the inventory-based estimates, supporting the validity of the modelling approach.

The simulated coarse woody litter production could not be directly compared with inventory-based estimates due to methodological issues. The litter and soil carbon model Yasso15 does not separate the deadwood, litter and soil organic carbon (SOC) pools (Tuomi et al.2011b; Didion et al. 2016). Based on earlier tests of the model, it can estimate the changes in the carbon pool of deadwood reliably given that input estimates are available (Didion et al. 2014; Hernández et al. 2017; Ziche et al. 2019). The simulated mean annual coarse woody litter input of the current study, 0.23-0.25 t C?ha?1?year?1, was in the same order of magnitude with other modelled estimates, 0.5 t C?ha?1?year?1from G4M and 0.08 t C?ha?1?year?1from EFISCEN, respectively (Repo et al. 2015). It is worth noting the definitions of coarse wood litter differ in different studies. The estimate of the coarse woody litter production does not compare directly to the actual coarse woody debris in the study area because the former is an input to deadwood pool and the latter is the pool. Nevertheless, the results demonstrate the impacts of residue extraction on the formation of coarse woody debris and therefore on the deadwood abundance in the long term. All in all, the simulated estimates of ecosystem services in the beginning of the simulation period were comparable with measurements, supporting the validity of the mapping framework.

The quality of the initial data on forest characteristics is also an important aspect of the performance of the mapping framework. A broad spatial coverage and comprehensive information on forest characteristics are the strengths of the MS-NFI data compared with coarser land use and land cover maps (Kangas et al. 2018). However, MS-NFI is more accurate on medium and large spatial scales rather than on individual grid cells. This is because the k Nearest Neighbour method averages stand volumes and site fertility classes levelling off extremes(Katila 2006; Haakana 2017). Furthermore, errors in the MS-NFI data are spatially autocorrelated (Katila and Tomppo 2001). Considering these limitations, the mapping framework is more suitable for assessing ecosystem services at landscape level rather than on individual forest management units.

The simulated stand age structure was a key determinant of the landscape level projections of ecosystem services in this study. The differences between the simulated and measured age class distribution reflected mainly the deviation of the MS-NFI data from the field plot measurements. The age class distribution was initially biased by overestimating the proportion of intermediate-aged stands because of the limitations of the k Nearest Neighbours method (Katila 2006; Haakana 2017). In the later simulation years, the forest management scenarios started to dominate over the initial status. In the current study, regular harvesting reduced the amount of fast-growing intermediate-aged stands risking carbon sequestration temporarily. The climate impacts of forest management depend, however, also on the life cycle of the wood products which was outside the scope of this study. In the future, the approach could be developed by setting a limit of harvest removal that could not be exceeded. Another option would be to implement knowledge of the actual harvest ages to the landscape simulation. The mapping framework could also be refined with multiple model runs to cover the potential variation in harvest regimes.

Implications for ecosystem service assessments

The mapping framework presented in this study contributes to bridging the gap between mapping and simulation modelling in the ecosystem service assessments of boreal forests. It incorporated new features in comparison with some existing tools for ecosystem service assessment (Nelson et al. 2009; Maes et al. 2012). Firstly,the couplings of biomass, litter and soil carbon cycles were accounted for by the modelling approach, like in many process-based models (Morales et al. 2005). As a result, the dynamics of forest carbon cycle were described more accurately than in tools utilizing simple,land cover -based proxies (Eigenbrod et al. 2010). Secondly, the simple structure of the mapping framework is an advantage compared with some detailed, computationally intensive forest simulators (e.g. Redsven et al.2004; Schelhaas et al. 2007; Rasinm?ki et al. 2009;) or process-based models (e.g. Bayer et al. 2015; Gutsch et al. 2018; Holmberg et al. 2019). The modular structure of the mapping framework enables its flexible development with new data and models in the future.Thirdly, the presented framework featured a stochastic development of forest age structure across the landscape,reflecting the variability of management regimes. This is a refinement in comparison with some decision-support systems applying fixed age classes (Frank et al. 2015).However, it is not possible to draw conclusions about the supremacy of the stochastic approach for age structure based on the current results.

In the future, indicators of water quality regulation(Wade et al. 2002; Huttunen et al. 2016) or functional diversity (Vihervaara et al. 2017) could be integrated to the framework to study their relationships with carbon cycle. Fine scale maps allow for investigating the spatial trade-offs and synergies between carbon sequestration,other ecosystem services and biodiversity. Fine scale layers of carbon stocks and fluxes could support sustainable land use planning when integrated into a spatial prioritisation system (e.g. Mikkonen and Moilanen 2013;Kukkala and Moilanen 2017). The large variety of forest types and management systems in the boreal zone sets a challenge for applying the mapping framework at broad spatial scales. For example, growth and yield models for old-growth and uneven-aged forests (Pukkala 2016), as well as litter and soil carbon models for organic soils(Ojanen et al. 2014), are few and require more development. It is also noteworthy that scenario-based applications alone are inadequate tools for finding optimal solutions for land use (M?nkk?nen et al. 2014). It would require simultaneous analysis of the alternative management regimes with multi-objective optimization tools(Eyvindson et al. 2018).

Conclusions

The mapping framework developed in this study integrated simulation modelling and spatially explicit, extensive data on forest characteristics. The approach produced reasonable estimates of the effects of bioenergy production on ecosystem services related to carbon cycle. It was suitable for assessing ecosystem services at the landscape level. Trade-off situations were observed between carbon sinks, wood production and coarse woody litter production as a result of continued harvest residue extraction for bioenergy production. The results demonstrated that stand age class distribution was a key driver of the simulated ecosystem service indicators across the study area. The framework contributes to bridging the gap between ecosystem service mapping and detailed simulation modelling in boreal forests. It allows for visualizing carbon stocks and fluxes as fine resolution maps to study their relationship with other ecosystem services and biodiversity. Fine scale maps of the impacts of forest management on carbon cycle could support sustainable land use planning. Future development of the framework includes integrating more detailed models and a wider variety of ecosystem service and biodiversity indicators to it.

Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s40663-021-00283-2.

Additional file 1:Figure S1 The development of the forest age class distribution in the Vanajavesi catchment area over the simulation period 2012-2100. Figure S2 The simulated (a) and measured(b) mean age class distribution in the study area in the beginning of the simulation period. The simulated and measured estimates represent the Vanajavesi catchment in 2012-2016 and the surrounding Kanta-H?me region in 2009-2013 (Natural Resources Institute Finland 2020),respectively.

Abbreviations

GPP: Gross Primary Production; MS-NFI: Multi-Source National Forest Inventory; NEP:Net Ecosystem Production; NBP: Net Biome Production;NPP: Net Primary Production;Rh:Heterotrophic Respiration; SOC: Soil Organic Carbon.

Acknowledgements

The authors wish to thank Prof. Miska Luoto (University of Helsinki),Dr. Mikko Savolahti and Dr. Niko Karvosenoja (Finnish Environment Institute), as well as the anonymous reviewers for their constructive comments on the manuscript.

Authors’contributions

AA, AR and JL outlined the original idea and study design. AA conducted the stand level forest simulations and AMA developed the stochastic landscape simulation. AA was a major contributor in interpreting the results and writing the manuscript with the support of AR. All authors read and approved the final manuscript.

Funding

AA has been supported by Maj and Tor Nessling Foundation through the grant “Coupling carbon sequestration of forests and croplands with ecosystem service assessments”(decision No. 201700251),LIFE+ financial instrument of the European Union (LIFE12 ENV/FI/000409, MONIMET) and the Academy of Finland Strategic Research Council project (SRC 2017/312559 IBC-CARBON).AR has been supported by the Academy of Finland through the grant“Trade-offs and synergies in land-based climate change mitigation and biodiversity conservation” (decision No. 322066).We declare that the funding bodies were not involved in the design of the study, collection, analysis or interpretation of data or writing the manuscript.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1Finnish Environment Institute SYKE, Biodiversity Centre, Latokartanonkaari 11, FI-00790 Helsinki, Finland.2University of Helsinki, Department of Geosciences and Geography, P.O. Box 64, FI-00014 Helsinki,Finland.3Natural Resources Institute Finland (Luke), P.O. Box 2,FI-00791 Helsinki, Finland.4VTT Technical Research Centre of Finland LTD, P.O. Box 1000, FI-02044 VTT Espoo,Finland.5Finnish Meteorological Institute, Climate System Research Unit, P.O.Box 503, FI-00101 Helsinki, Finland.

Received: 9 April 2020 Accepted: 12 January 2021