Sources of uncertainty

Sources of uncertainty are identified adjusting the uncertainty taxonomy described, e.g. in Kujala et al. (2013) and McGlynn et al. (2022). Uncertainty elements in the simulations were: (1) inputs: (1a) forest initial state and soil fertility and (1b) GCM uncertainty representing structural uncertainty in the used GCMs; and (2) parameters: (2a) model parameters and (2b) harvest target level. The simulation settings are shown in Fig. S2. The input uncertainties varied spatially in each region, while the parametric uncertainties varied only between simulations, i.e. they remained constant over different regions and through simulation time periods. Although the climatic conditions vary spatially, they remain constant under individual GCMs; thus, they were categorised as input uncertainty in this study. The drivers that were set as constant in the simulations were the climate and harvest scenarios, which were categorised as human decision uncertainty. See Table 3 for different types of uncertainty elements.

Uncertainty elements in the simulation

Epistemic uncertainty is due to the lack of knowledge of the system in respect to quantities and processes within the system. Aleatory uncertainty arises because of the unpredictable, random nature of the physical system under study. Human decision uncertainty arises from subjective human preferences and beliefs

In the spatially explicit, high-resolution direct Monte Carlo simulations, computational effort is heavy even when using the segmented input variables as computation units: The simulations are repeated multiple times in each scenario for each computational unit. With the 28 million segments covering the forest area of this study, the computational effort needed to be reduced.

The MS-NFI pixel level results are outputs of an improved k nearest neighbours (ik-NN) model, which is based on the NFI field measurements and spatially projected satellite data (Mäkisara et al. 2019). These pixel level variable values can be considered as realisations of the real variable values plus random modelling error. The segment level data used as the initial state variables included also the random segmentation errors and served here as the erroneous initial state variable value population.

In the region level analysis, the focus is on the region level total and average output results, not on the pixel or segment level. For such analysis, each initial value set of the simulations has to be a representative sample of the population of the initial values of that region. To reduce the computational effort in this study, region level simulations were performed for a random set of 20 000 pixels in the region. The pixels were sampled with replacement from the segmented data population using the proportion of the distinct segment areas from the total study area as weights in order to generate a representative sample of the pixels.

Measured error estimates were not available for the mean age of trees in a pixel. Thus, to produce a realistic approximation of the age data precision, a rough estimate for the age uncertainty based on the experience of the authors was used instead. The mean age of trees in the pixel was assumed to be known more precisely for the younger trees than the older trees. For each pixel, they were sampled using a normal distribution with age-dependent standard deviation: xage,j∼Nμage,j,σage,j2, where μage,j is the MS-NFI mean tree age of the pixel j and σage,j=0.1μage,j. Thus, the 95% probability range of 10 years old trees is approximately 8–12 years, while for 100 years old trees it is 80–120 years.

Pixel level site type is also an estimate in the initial state data. Uncertainty in the site type was simulated by re-sampling site type of each sampled pixel at the beginning of each simulation. The model for the site type probability distributions followed the model given in Haakana et al. (2022). The pixel’s site type probability distribution depends on the MS-NFI based site type and sampled mean tree height, mean tree breast height diameter, basal area and proportion of pine trees. The probability distribution was estimated for each pixel according to its structural variable values. The resulting probability distribution was used to sample a new site type for that pixel. The sampled site type was kept constant over the whole simulation period.

The initial values of PREBAS simulations did not include information about the initial state of the soil carbon stock on mineral soils. The soil carbon processes were estimated using the YASSO model combined with the litter outputs from the PREBAS model. In the absence of the measured data of the initial state, it was modelled with assumption of a steady state. Here, the steady state was estimated for each simulation i separately, starting from the random initial state variables described above. PREBAS model with randomly ordered repetitions of realised, historical harvest levels and real local weather data from the years 2015–2021 was run until steady state was reached. The steady state, thus, depends on the sampled initial values and model parameters. The same initial state was used for projections of different harvest and climate scenarios.

Model parameters include the parameters used in different sub-model components attached in the PREBAS model: CROBAS, PRELES and YASSO. Samples of these parameters have been estimated and validated with measured data in previous studies (Minunno et al. 2016, 2019; Viskari et al. 2021). These samples of parameter values were used as the parameter populations, from which random sample sets for each simulation were re-sampled with replacement for each simulation i.

In the PREBAS model, the initial value of the species-dependent volume is a function of the forest structural values given in the MS-NFI data and estimated crown height. Crown height is not included in the MS-NFI data; thus, it has been estimated using empirical equations (Sharma et al. 2017). However, uncertainty of the crown height estimate was not available for this study; thus, a rough estimate for the crown height uncertainty based on the experience of the authors was used instead. Here, the crown height uncertainty was simulated by sampling the crown height factor from normal distribution, ccrownheight,i∼N(1,0.12) and using it to multiply the estimated crown height. Here the crown height was assumed to vary between 80 and 120% of the estimated crown height with 95% probability.

The PREBAS and YASSO models are based on assumption of forests growing on mineral soil. However, also forests located on drained organic soils were included in the simulations. The segments located on drained organic soils have been classified according to the site type to correspond with mineral sites of similar fertility, and the organic soil impact on growth was modelled based on this classification. Uncertainty of the classification of mineral and drained organic soils was not simulated in this study.

Uncertainty in the average organic soil emissions was simulated by sampling emission coefficients, EFet∼N(μet,σet2), for which the mean and standard deviations for the emission type (et) are given in Table 1. The sampled emission coefficients remained constant during the simulation period and over all sampled pixels and regions of simulation i. The soil carbon stock of forests on organic soils was sampled from xsoilC,peatland∼N(543400,185002) [kg C ha-1] for all site types (Turunen and Valpola 2020).

In the climate scenarios RCP2.6, RCP4.5 and RCP8.5, the weather data uncertainty was a result of the variation within different GCM’s: CanESM2, CNRM, GFDL, HadGEM2 and MIROC. In each simulation i, one of the five climate models was chosen randomly with equal probability. The same GCM was used over all the regions. This approach allowed computationally efficient validation of the effect of GCM based uncertainty within a limited number of simulations.

Historical data of the whole country level harvests were estimated based on information given by Natural Resources Institute Finland (2021). The uncertainties in round wood and energy wood statistical harvest levels were simulated by sampling the target level using normal distribution with mean value given in statistics or in harvest scenario specific projected target level, and standard deviation as 2% of the mean (Peltoniemi et al. 2006). The sampled whole country harvest levels were allocated to region level according to the region harvest level proportion in historical data.

Separate uncertainty source components given in this section were sampled independent of each other; thus, the possible correlations between different components of uncertainty sources were ignored. In the input data uncertainty, the possible spatial correlation of site type index and age were ignored, as the homogeneous segments were considered spatially independent. Also, the PREBAS model relies on the assumption of no interaction between segments. Classification of forest land class (forest land and poorly productive forest land) in MS-NFI data and classification of forests location to mineral soil or to different types of organic soils were based on high resolution spatial data and were assumed accurate in this study.

The most significant uncertainties that were not included in this study, are forest disturbances such as wind damages, forest fires, snow damages and biotic risks. They are not included in the current PREBAS version, and overall, spatially explicit projections of such events are difficult to obtain.