Materials
Study area
The study area (Fig. 1) lies in Krødsherad municipality, southeastern Norway. It is a typical boreal forest dominated by Norway spruce (Picea abies (L.) Karst.), Scots pine (Pinus sylvestris L.), and, to a less extent, birch species (Betula pendula Roth and Betula pubescens Ehrh). The area consists of 3324 managed forest stands spanning approximately 50 km2.
Field data
A total of 116 circular area plots (232.9 m2) were distributed systematically within three strata: young forest (39 plots), mature forest with poor site quality (38 plots, site index < = 11), and mature forest with good site quality (39 plots, site index > 11). The plots were measured in two campaigns, approximately 15 years apart. The first campaign was conducted in 2001 followed by a second when all plots were revisited in 2016 and 2017. For the second campaign that spanned two years, 2016 was set as reference year. The diameters of all trees above 10 cm (4 cm in the young forest stratum) were measured with a caliper, and the heights of a relascope-based subsample of around 10 trees per plot were measured with a hypsometer. The following variables were calculated: mean diameter (\(D\)), mean height (\(H\)), dominant height (\({H}_{dom}\)), and number of trees per hectare (\(N\)). Using allometric models [41], total living biomass per hectare (\(BMS\)) was estimated.
Airborne laser scanner data
The ALS acquisition temporally overlapped the field data collection, the entire study area having been scanned under leaf-on conditions in 2001 and 2016. The field surveys and ALS acquisitions are described in greater detail in [42] and [43].
The ALS point clouds were normalized and tessellated with a grid of 15.26 m\(\times\)15.26 m cells having the same area as the field plots. For each cell, metrics describing the vertical distribution of the laser returns were computed:
- Height percentiles \({H}_{ALS\_10}\), \({H}_{ALS\_20}\), …, \({H}_{ALS\_90}\) corresponding to the 10th, 20th, …, 90th percentiles
-
Mean height (\({H}_{ALS\_mean}\))
-
Cumulated densities \({D}_{ALS\_1}\), \({D}_{ALS\_2}\), …, \({D}_{ALS\_10}\) as proportions of points above ten equally spaced height levels from 1.3 m to the 95th height percentile.
The set of ALS metrics were calculated for the 116 georeferenced sample plots as well.
Forest stand map
A forest stand map was obtained from the forest management inventory carried out in the area in 2018. Each stand had the following attributes: site index (\(SI\)), stand age (\(AGE\)), and species proportions by volume recorded, which determined the dominant species (\(SP\)). The stand attributes were obtained using a combination of projections of the old stand map and updates using photointerpretation of aerial imagery.
Satellite imagery
Landsat satellite imagery were used to detect large disturbances, such as harvest and thinnings, at stand level. We used the temporal segmentation algorithm LandTrendr [44], implemented with Google Earth Engine [45]. LandTrendr has been widely used and shows good performance in similar environments [46]. A disturbance map was created based on the ALS tessellation and recording the year of disturbance.
Tree allometry dataset
Marklund’s [41] allometric models have been extensively employed in Norway and Sweden for many decades. They are used to predict biomass for individual tree components. The original publication reported the estimated model parameters but did not include any materials on the errors associated with the parameters. Using the original set of observations, we refitted the models and estimated the covariance matrices for the parameters, which are needed for the statistical inference.
The dataset consisted of 1281 individual trees (546 spruce, 494 pine, and 241 birch) with measured heights and diameters. The response variables were biomass of stem wood (\(SW\)), branches (\(BR\)), dead branches (\(DB\)), bark (\(SB\)), stump (\(SU\)), foliage (\(FL\)), fine roots (\(RF\)), and coarse roots (\(RC\)). The number of observations for each biomass component varies as shown in Table 1.
Table 1
Number of observations for each biomass component by species
|
\(SW\)
|
\(BR\)
|
\(DB\)
|
\(SB\)
|
\(SU\)
|
\(FL\)
|
\(RF\)
|
\(RC\)
|
Spruce
|
521
|
540
|
525
|
521
|
323
|
540
|
324
|
333
|
Pine
|
471
|
486
|
472
|
471
|
305
|
485
|
305
|
314
|
Birch
|
218
|
235
|
221
|
218
|
0
|
0
|
0
|
0
|
Climate data
The climate data for this study have been acquired using the Frost API of the Norwegian Meteorological Institute [47]. The data are licensed under Norwegian license for public data (NLOD) and Creative Commons 4.0 BY. Historical weather data were retrieved by identifying the nearest weather station and retrieving timeseries of daily precipitation and temperature data. For each year from 2001 to 2016 the following climate variables were calculated: annual precipitations (mm), mean annual temperature (°C) and the mean difference between maximum and minimum monthly temperatures (°C). Another set of climate variables was calculated with data starting from 1957, the earliest recorded year to 2000. For these older observations, the three variables were averaged across the 43 years.
Methods
Overview
We start with a brief overview of the proposed methodology. First, we introduce the five variables of interest to be estimated at the stand level: \({\varDelta C}_{AGB}\), \({\varDelta C}_{BGB}\), \({\varDelta C}_{litter}\), \({\varDelta C}_{deadwood}\), and \({\varDelta C}_{SOC}\). They represent the change in carbon mass expressed in Mg ha− 1 within the following pools: above-ground biomass (AGB), below-ground biomass (BGB), litter, deadwood, and soil organic (SO) carbon.
We use an indirect method to estimate change [31], by taking the difference between estimates of C stocks at the start and end of the time period (i.e., \({\varDelta C}_{AGB}={C}_{AGB|2016}-{C}_{AGB|2001}\)). Moreover, in lack of ground observations we rely on a model-based approach where the inference is based on the properties of the models involved in the estimation [48].
The carbon stocks in the living biomass pools (\({C}_{AGB}\) and \({C}_{BGB}\)) were estimated for both points in time using area-based ALS models. The carbon stocks in the dead biomass pools (\({C}_{deadwood}\) and \({C}_{litter}\)) and \({C}_{SO}\) pools, however, are accumulations sourced from the biomass pools that go through a decomposition process. To estimate their levels, we need to: 1) approximate the litter and deadwood production, and 2) simulate the decomposition process. For the first part, we calculated the litter production based on the living biomass stocks and active silvicultural treatments. We considered three mechanisms that generate litter: annual litter turnover, annual mortality, and excess litter resulting from harvest. The decomposition process was simulated using the Yasso15 soil model [49].
The estimation processes involve several models (see section 2.2.2) that are linked together in a chain of predictions (see section 2.2.3), the outcome being carbon stocks in the five pools at cell level. For the soil-related pools (\({C}_{deadwood}\), \({C}_{litter}\), and \({C}_{SO}\)) the carbon stocks are calculated year by year starting from 2001, through 2016, each time carrying through the accumulated carbon from the previous year and integrating new yearly litter.
The stand level estimates of carbon change in each pool are obtained by averaging predictions over the cells within each stand. To estimate the uncertainty, we used a Monte Carlo approach, by sampling repeatedly from the models’ parameters distributions using their estimated means and covariance matrices and assuming joint normal distributions (section 2.2.3.1).
Models
Yasso15 model
The Yasso15 model [49] partitions soil in five chemical compartments. Four of these compartments belong to the decomposing litter: celluloses (\(A\)), sugars (\(W\)), wax-like compounds (\(E\)), and lignin-like compounds (\(N\)). The fifth compartment is humus (\(H\)) as the end of the decomposition process. The carbon accumulated in humus makes up the soil organic (SO) carbon pool. The Yasso15 model is formulated in terms of decomposition rates (Fig. 2). Each litter compartment decomposes to either: another litter compartment, humus, or is emitted as CO2. The carbon in the humus compartment is only emitted as CO2. The decomposition rates depend on three climate variables: annual precipitations (mm), mean annual temperature (°C) and the mean difference between maximum and minimum monthly temperatures (°C).
Since the flow rates of the Yasso15 model depend only on the climate variables, it allows to separate carbon from different sources and trace their decomposition independently. We used this property to separate the decomposition of the initial carbon stocks (e.g., stocks at the beginning of the timeframe) from the new carbon being stored from the litter produced during the timeframe of interest. Also, we treated separately carbon that originates from normal turnover and harvest (\({C}_{litter}\)) from carbon sourced in mortality (\({C}_{deadwood}\)).
Allometric models
The allometric models for tree component biomass are based on Marklund [41]. The models were refitted using the original set of field observations and the parameter covariance matrix was estimated together with the mean values. The models are of three forms:
(1)\(\text{ln}\left(comp\right)={\beta }_{0}+{\beta }_{1}\frac{D}{D+k}\)
(2)\(\text{ln}\left(comp\right)={\beta }_{0}+{\beta }_{1}\frac{D}{D+k}+{\beta }_{2}\text{l}\text{n}\left(H\right)\)
(3)\(\text{ln}\left(comp\right)={\beta }_{0}+{\beta }_{1}\frac{D}{D+k}+{\beta }_{2}H+{\beta }_{3}\text{ln}\left(H\right)\)
where \(comp\) is one of: \(SW\), \(BR\), \(DB\), \(SB\), \(SU\), \(FL\), \(RF\), \(RC\).
Table 2 shows which model form was used for each component and species.
Table 2
Summary of the allometric models for biomass components
|
\({S}{W}\)
|
\({B}{R}(+{F}{L})\)
|
\({D}{B}\)
|
\({S}{B}\)
|
\({S}{U}\)
|
\({F}{L}\)
|
\({R}{F}\)
|
\({R}{C}\)
|
Spruce
|
Eq. 3, k = 14
|
Eq. 3, k = 13
|
Eq. 3, k = 18
|
Eq. 3, k = 15
|
Eq. 1, k = 17
|
Eq. 2, k = 12
|
Eq. 1, k = 12
|
Eq. 1, k = 8
|
Pine
|
Eq. 3, k = 14
|
Eq. 2, k = 10
|
Eq. 3, k = 10
|
Eq. 2, k = 16
|
Eq. 1, k = 15
|
Eq. 3, k = 7
|
Eq. 1, k = 10
|
Eq. 1, k = 9
|
Birch
|
Eq. 2, k = 11
|
Eq. 1, k = 10
|
Eq. 3, k = 30
|
Eq. 2, k = 14
|
Use Pine
|
\(\frac{0.011}{0.52}\times SW\)
|
\(\frac{0.042}{0.52}\times SW\)
|
\(\frac{0.042}{0.52}\times SW\)
|
The k values are used in Eq. 1–3
|
Area based models
The area-based models were fitted using the 116 field plots for which both field-calculated biophysical variables and ALS metrics were available. All area-based models were selected using the Bayesian information criterion, restricting the maximum number of parameters to five and the maximum variance inflation factor to five. The models were time-invariant [30], fitted on observations from both points in time, and using a dummy variable: \(T=0\) for 2001, and \(T=1\) for 2016. The models had the following form:
Mean diameter model:\(\text{ln}\left(D\right)={\beta }_{0}+{\beta }_{1}{D}_{ALS\_4}+{\beta }_{2}{H}_{ALS\_mean}+{\beta }_{3}T\)
Dominant height model:\({H}_{dom}={\beta }_{0}+{\beta }_{1}{D}_{ALS\_2}+{\beta }_{2}{D}_{ALS\_8}+{\beta }_{3}{H}_{ALS\_80}+{\beta }_{4}T\)
Biomass model:\(\text{l}\text{n}\left(BMS\right)={\beta }_{0}+{\beta }_{1}{D}_{ALS\_2}+{\beta }_{2}{H}_{ALS\_80}+{\beta }_{3}T\)
Based on the same dataset a simple model was fitted to predict mean tree height (\(H\)) using \({H}_{dom}\):
Mean height model:\(H={\beta }_{0}+{\beta }_{1}{H}_{dom}+{\beta }_{2}T\)
We needed to predict both mean and dominant height since the growth models work with dominant heights and the biomass component models are for individual trees, thus using the mean height.
Models with unaccounted uncertainty
There are several external models for which the authors published only the mean parameter estimates, and thus we could not account for the uncertainty in the estimated parameters.
The diameter growth models published by Blingsmo [50] were fitted on plot data from the Norwegian national forest inventory (NFI). The growth period was on average five years, and the number of observed periods for each species was: 1385 for spruce, 1292 for pine, and 662 for birch. The diameter growth models were fitted separately for each species with the following dependent variables: diameter (\(D\)), site index (\(SI\)), number of trees per hectare (\(N\)), dominant height (\({H}_{dom}\)), and age (\(AGE\)).
The dominant height growth models were published by Sharma, Brunner [51] for spruce and pine and Eriksson, Johansson [52] for birch. The dependent variables were \(SI\) and \(AGE\).
The litter turnover rates and the AWEN partition of tree biomass components were used as fixed values, with unaccounted uncertainty. The carbon content of biomass was also a constant, i.e., 50%.
Estimation process
Computing the estimated change in soil and biomass carbon involved combining several models including ALS area-based models, allometric models, growth models and the Yasso15 soil model which were described in the previous section. For a better overview, we grouped the models in several functional blocks (Fig. 3). The core process is illustrated in Fig. 4 and Fig. 5 shows how the final predictions for each of the five pools are calculated.
The mean tree in each cell is characterized by its diameter (\(D\)), height (\(H\)), and species and it was predicted using ALS area-based models (Fig. 3. A). \(H\) was predicted in two steps, with the dominant height model as intermediary. The growth models are illustrated in Fig. 3. B. Again, \(H\) is obtained via \({H}_{dom}\).
Finally, Fig. 3. (C) illustrates the litter calculation. Starting with the mean tree, the biomass components are predicted using the Marklund models. Next, depending on the litter source and tree species, a certain fraction of each component is retained as litter. We consider three litter sources: normal turnover, mortality, and harvest. The biomass fractions for the normal turnover are shown in Appendix A. For mortality, the fractions are 100% of all components, and in case of harvest, 95% of the stem wood is extracted (i.e., 5% left as litter), while the rest of biomass components are left in the forest and turn to litter entirely (100%). The yearly mortality rate was determined by the difference in the estimated number of stems (see Fig. 5. (E)) at the two points in time, and in the case of non-decreasing stem number a constant rate of 0.4% was assumed yearly [53].
Finally, the litter originated from each biomass component was partitioned into four compartments according to the chemical composition: \(A\), \(W\), \(E\), \(N\) (see Appendix B).
The core prediction process is illustrated in Fig. 4. It consists in estimating the mean tree at key points in time, calculating the litter associated with that, obtaining the yearly litter production by interpolating the litter quantities for the years in between, and finally use the Yasso15 model in a chain of predictions.
Since evolution of the growing stocks depends on whether the forest grew undisturbed or there was a disturbance event, we considered two scenarios: undisturbed forest growth, and disturbance detected. The forest was assumed to be disturbed if the estimated biomass from the ALS survey (\(BMS\) – area based) decreased between 2001 and 2016. In this case, the disturbance map would provide the approximate year of the disturbance event. If no year was recorded, then the event was assumed to have happened in the middle of the time interval. In the simpler case of undisturbed growth, the yearly litter was calculated by interpolating linearly between the litter quantities associated with the mean tree at the start and the end of the time interval. If disturbance was detected, two additional mean trees were predicted: the mean tree right before the event, and the mean tree right after the event. The “before” tree was predicted using the growth models with the tree in 2001, and the “after” tree was predicted using inverse growth models with the tree in 2016. As the growth models are difficult to invert analytically, an approximation search algorithm was used, where the search interval was recursively halved. The tolerance for both height and diameter predictions were set to be less than 10− 3 m. In this scenario we have four points in time with determined growing stocks size, and corresponding litter production. We interpolate linearly for the years in between 2001 and the year of the event, and from the following year to the end in 2016. For the year of the disturbance event, we assume an active silvicultural treatment (thinning or harvest), and the normal litter production is supplemented by the vegetal residuals typically left on site (i.e., all but 95% of the stem wood).
The initial soil carbon quantities (at the start of the timeframe; 2001) denoted here as “old” soil carbon were approximated using simulations in two stages. First, the long-term soil carbon was calculated by running the model iteratively with a fixed litter input until the carbon quantities in the chemical compartments reach an equilibrium. The litter input was determined using average values of the growing stocks in 2001 for the entire study area separately for each tree species and site index (Table 3). In absence of empirical soil carbon observations or knowledge of the old history of the stand, this ensured that all stands with similar productivity have a common soil carbon baseline.
Table 3
Long-term soil carbon values (Mg ha− 1) by species (SP) and site index (SI)
\(SP\)\(SI\)
|
\(SP\)
|
\(SI\)
|
\(A\)
|
\(W\)
|
\(E\)
|
\(N\)
|
\(H\)
|
Spruce
|
6
|
6.97
|
0.73
|
0.68
|
17.74
|
35.09
|
8
|
5.45
|
0.57
|
0.53
|
13.88
|
27.43
|
11
|
7.10
|
0.75
|
0.69
|
18.07
|
35.77
|
14
|
7.16
|
0.76
|
0.70
|
18.23
|
36.06
|
17
|
7.94
|
0.84
|
0.78
|
20.21
|
39.99
|
20
|
8.17
|
0.86
|
0.80
|
20.78
|
41.12
|
23
|
8.36
|
0.88
|
0.82
|
21.29
|
42.11
|
Pine
|
6
|
4.56
|
0.48
|
0.51
|
11.33
|
23.02
|
8
|
5.64
|
0.60
|
0.63
|
13.98
|
28.43
|
11
|
7.35
|
0.78
|
0.81
|
18.21
|
37.04
|
14
|
8.35
|
0.88
|
0.92
|
20.68
|
42.05
|
17
|
9.98
|
1.05
|
1.09
|
24.73
|
50.28
|
20
|
11.81
|
1.25
|
1.29
|
29.25
|
59.45
|
Birch
|
8
|
5.15
|
0.56
|
0.62
|
13.54
|
26.42
|
11
|
7.34
|
0.79
|
0.90
|
19.33
|
37.66
|
14
|
7.50
|
0.81
|
0.91
|
19.74
|
38.49
|
17
|
8.67
|
0.94
|
1.06
|
22.84
|
44.49
|
20
|
9.41
|
1.02
|
1.15
|
24.80
|
48.31
|
23
|
14.06
|
1.52
|
1.74
|
37.16
|
72.27
|
In the second stage, the evolution of soil C was calculated for the current silvicultural cycle which was assumed to have started with a clear cut and the soil carbon values in Table 3. The stand \(AGE\) in 2001 determined how many years have passed in the current silvicultural cycle. The soil carbon model was applied for each year in the current cycle until 2001, with the values for the yearly litter input being linearly interpolated between 0 and the ones calculated for the year 2001.
We expect the soil carbon initialization to be a coarse approximation of the soil carbons stocks in the year 2001, so we kept this value separated from the predicted accumulation during the timeframe that is based on several good quality data sources and models. This means that: (1) we traced the old soil carbon separately through the 15-year timeframe, with no yearly litter input, and (2) the soil carbon starts with 0 initial values for the within timeframe process (Fig. 4). Finally, the old carbon stocks in the year 2016 may be added to the new carbon accumulated in the timeframe, the result being identical to having the timeframe processing initialized with the old carbon in 2001.
The final steps to obtain carbon stocks in the five pools are shown in Fig. 5 (E – for AGB, and BGB; F – for litter, deadwood, and SO). The \({C}_{AGB}\) and \({C}_{BGB}\) stocks are calculated using the related mean tree biomass components (i.e., \(RF\) and \(RC\) belong to BGB and the rest to AGB), and then scaling up to per hectare values using the total \(BMS\) predicted with the area-based model. The number of trees per hectare is calculated by dividing \(BMS\) to the biomass of the mean tree. \({C}_{litter }\) and \({C}_{deadwood}\) are calculated separately using the process in Fig. 4. \({C}_{litter }\)accounts for the normal litter turnover plus the harvest residues, and \({C}_{deadwood}\) accumulates carbon sourced in mortality. \({C}_{litter }\) and \({C}_{deadwood}\) are calculated by summing up the carbon in the litter chemical compartments: \(A+W+E+N\). Finally, \({C}_{SO }\)consists of the humus (\(H\)) compartments of both \({C}_{litter }\) and \({C}_{deadwood}\). Table 4 shows the connection between the Yasso15 chemical compartments, and the soil carbon pools.
Table 4
Summary of the Yasso15 chemical compartments, carbon input origin, and the soil carbon pools
|
Timeframe
|
Historical
|
|
To. + Hv.
|
Mort.
|
To.
|
Mort.
|
A
|
\({C}_{litter}\)
|
\({C}_{deadwood}\)
|
\({C}_{litter\left(old\right)}\)
|
W
|
E
|
N
|
H
|
\({C}_{SO}\)
|
\({C}_{SO\left(old\right)}\)
|
The abbreviations are turnover (To.), harvest (Hv.), and mortality (Mort.)
|
Model based estimation and uncertainty
As shown in the previous sections, predicting the soil carbon change in a cell involved the iterative use of the Yasso15 soil model and a series of interconnected models to calculate the yearly litter. The resulting soil carbon change predictions were then aggregated at the stand level. Since tracing the errors analytically would be extremely tedious, we used parametric bootstrapping [54]. This is a Monte Carlo-style method where the parameter values of each model were iteratively sampled from their estimated distributions, each time recalculating the whole chain of predictions to a new outcome. For the five models that we fitted ourselves (i.e., area-based biomass, dominant height, mean height, diameter, and biomass components - Marklund) we sampled from the joint parameter distributions defined by estimated means and variance-covariance matrices [55]. For Yasso15, the Finish Meteorological Institute provided us with a readily generated sample of 10000 pairs of parameter values. For each parameter sample the entire sequence of predictions was recalculated to a different outcome (i.e., carbon change in the five pools). The errors were thus approximated by the sampling distribution of the change estimators.
To assess the contribution of an individual model to the total error, we ran separate Monte Carlo simulations for each of the six models, where parametric bootstrapping was applied to the parameters of one model at a time, keeping the parameter of the rest of the models fixed at their estimated mean. This type of analysis does not account for the interaction between the models, the separate error components do not add up to the errors calculated for all models at once, so we report their relative size as percentage of the total error.