Estimating Carbon Storage in Tree Stem Biomass of Mixed Boreal Forest Using UAV LiDAR Data

1. Introduction

Forests are the most important terrestrial ecosystems, providing accumulation and long-term storage of biogenic carbon in the form of timber and soil horizons rich in organic matter [1,2]. Inventorying forests' carbon stocks and assessing their spatiotemporal variability is a critical task in the context of climate change and the search for ways to adapt to it, as noted at both international [3] and national [4] levels. Since 2022, Russia has been implementing the project "National system for monitoring carbon pools and greenhouse gas fluxes in Russia" ("Rhythm of Carbon"), one of the main goals of which is the assessment of carbon stocks in forest ecosystems, which is already being carried out at the developing network of permanent monitoring stations [5].

The most significant carbon pool of boreal forests is the phytomass of stem wood, which accumulates up to 20–30% of carbon in such ecosystems [1,6]. Traditionally, the inventory of stem wood stocks (and consequently carbon) in forests is carried out using field forest inventory methods on standard sample plots of 50×50 to 100×100 meters. Such data requires further scaling in space and time to inventory of carbon stocks at local and regional levels [7,8,9].

Upscaling information obtained on the ground is conveniently done using remote sensing satellite data: information on location, condition, and predominant tree species is available in the visible and multispectral ranges [7,10,11,12,13,14]; in the lidar and radar ranges, information on tree height and crown sizes is available [15,16,17,18,19]. However, the resolution of open satellite data (Landsat, Sentinel) is often insufficient for accurate interpretation of ground information both spatially and temporally, especially at the level of individual forest inventory plots of small size, and commercial imagery is expensive and does not always sufficiently cover the territory for which upscaling is carried out [20,21].

The link between ground and satellite data in the process of their spatial and temporal interpretation can be effectively achieved using unmanned aerial vehicles (UAVs). They allow for the identification of small-scale (at the level of individual trees) areas of forest cover degradation and recovery, obtaining information on forest productivity, predominant tree species, phenological stage, fires, and pest spread. In addition, LiDAR sensors installed on UAVs enable the assessment of the vertical structure of forest communities, morphological characteristics of the stand, features of the underlying terrain, stratification, tree growth dynamics, and the condition of individual trees. It should be emphasized that data at the level of individual trees (sub-centimeter resolution) can only be obtained using UAVs, which significantly expands the capabilities of ground-based forest inventory work for upscaling, including with the help of satellite data [22,23,24].

Tree height, assessed using LiDAR UAVs, is a key parameter in evaluating woody biomass stocks; when combined with considerations of tree species, crown projective cover, and trunk diameter, it is sufficiently accurate for calculating carbon stocks in them [25,26], especially in the presence of allometric models with carefully identified parameters. Determination of tree heights is performed using a three-dimensional point cloud of laser reflections (LiDAR), obtained by recording the direction and travel time of the laser scanner beam from source to receiver, with a root mean square error of about 1–10 cm [27,28], which is comparable to or more accurate than ground-based assessment.

Remote LiDAR data require careful processing and ground verification: errors in estimating the number and height of trees often arise under conditions of rugged terrain, high canopy closure, low tree heights, and insufficient density of the point cloud of laser reflections or small overlap area of images for photogrammetric processing. In [29], it was noted that measuring tree height based on LiDAR surveying tends to underestimate the height of tall trees and overestimate the height of short trees; in [30], it was shown that some algorithms for processing LiDAR data tend to underestimate or overestimate tree heights; according to our earlier studies [31], when assessing the height of tree cover in an oligotrophic bog area (low-growing pines with an average height of up to 3 meters) using LiDAR surveying, a systematic underestimation of tree height was noted. Nevertheless, approaches based on the use of LiDAR surveying contribute to obtaining sufficiently accurate and reproducible results in the assessment of morphological properties of tree cover [27,28,30].

This work is dedicated to the verification of stem wood carbon stocks obtained from LiDAR UAV surveying by means of ground-based forest inventory in a mixed boreal forest site in Western Siberia.

2. Materials and Methods

2.1. Study Area

The research was conducted at the mixed forest site of the international field station "Mukhrino," located on the left bank of the Ob and Irtysh rivers, 20 km southwest of the city of Khanty-Mansiysk, Russia [32], in June 2023 (Figure 1). The mixed forest site, measuring 1.5×4.5 km (center coordinates: N60.897437 E68.711363), is a permanent plot for ground-based carbon stock assessment, which will be scaled up (together with other plots as part of the development of the observation network) to the territory of Russia within the "Rhythm of Carbon" project using satellite remote sensing data and mathematical modeling. The study area is located in the middle taiga of the boreal biogeographical zone, and the mixed herbaceous-green-moss forest with a canopy closure of 70-80% consists of two layers. In the first tree layer, aspen (Populus tremula), birch (Betula pubescens), siberian pine (Pinus sibirica), fir (Abies sibirica), and spruce (Picea obovata) predominate. The understory layer includes fir, siberian pine, and spruce. The undergrowth mainly comprises siberian pine, fir, spruce, and aspen. The total number of live undergrowth ranges from 7000 to 19000 plants per hectare. The stand is of III–IV site class, reaching heights of 20–25 m and trunk diameters of 30–50 cm. The understory is poor in species composition, mainly represented by rowan, and more rarely by rosehip. The number of understory plants ranges from 400 to 3800 per hectare. In the ground cover, taiga small grassy plants dominate – Linnaea borealis, Maianthemum bifolium, Gymnocarpium dryopteris, and boreal species of green mosses Hylocomium splendens and Pleurozium schreberi. Less abundantly, there is Vaccinium myrtillus, Equisetum sylvaticum and species of clubmosses like Diphaziastrum complanatum and Lycopodium annotinum. The forest represents a climax stage of development on well-drained sandy loam and loamy podzols. These mixed forests are typical for the middle taiga of Western Siberia and are often associated in the landscape with excessively moist types, such as swampy birch-sedge (Carex globularis) sphagnum (Sphagnum angustifolium) and horsetail-small grassy forests, as well as at the borders with bogs, the so-called "ryams", which often form continuous landscape gradients. Ryam is a type of forested bog characterized by a shrub layer of Ledum palustris and Chamaedaphne calyculata, with the participation of Vaccinium myrtillus, and a tree layer of 6–10 meters high Pinus sylvestris with a small admixture of Siberian cedar (P. sibirica), birch (Betula pubescens), and predominance in the ground cover of Sphagnum angustifolium and Sphagnum divinum.

2.1. Forest Taxation And Allometric Models

In the mixed forest, 5 sample sites were established (the coordinates of the site centers are presented in Table 1) for conducting forest inventory, which included: a complete enumeration of trees with separation by layers and represented species, measurement of tree heights (using a SUUNTO PM-5/1520 hypsometer) and diameters at breast height (at 1.3 m using a diameter caliper), assessment of above-ground phytomass, as well as evaluation of wood density in model trees considering layers and species.

The aboveground phytomass of the model trees (one average tree for height and diameter in each of the two layers of each sample plot) was determined for the overstory layer by direct weighing of the sawn trunk using KERN CH50K scales; and for the trees of the understory layer – by calculation (by allometric model).

To achieve this, the trunk volume was calculated (equation 1) and the wood density was measured (equation 2) in disks with a thickness of 3–4 cm, cut out from the middle of two-meter trunk sections (at odd meters).

where: V_g – trunk volume, m³; l – section length, m (all sections except the top had a length of 2 m); s_i – cross-sectional area in the middle of the two-meter section (at odd meters, the index i corresponds to the height in meters, for example, s₁ – area at a height of 1 m, s₃ – at a height of 3 m, etc.), m² (measured for the cut disks); s_top – cross-sectional area at the base of the tree top, m²; l_top – length of the tree top, m.

$$\mathsf{\rho }={\displaystyle \frac{m}{V_d}},$$

(2)

where: ρ – density of the cut disk, kg/m³; m – mass of the disk, kg; V_d – volume of the disk, m³. The mass of the disk was determined using KERN 440-47 scales. The volume of the disk was determined by four measured diameters (in different directions, averaged) and by eight measurements of thickness (also averaged). The mass of each two-meter section was obtained by multiplying the density of the corresponding disk ρ by the volume of the corresponding two-meter section; the sum of the masses of all sections was equal to the trunk phytomass.

The weighted average wood density for each of the five studied plots was calculated using equation 3.

$${\mathsf{\rho }}_w={\displaystyle \sum _{i=1}^n{\mathsf{\rho }}_i\cdot w_i},$$

(3)

где: ρ_w – is the average weighted density of wood in the area, kg/m³; ρ_i – is the density of the wood of the i-th species of trees in the area, kg/m³; w_i – is the dimensionless weight coefficient calculated as the ratio of the mass of the trunks of the i-th species of trees to the total mass of the trunks of all trees in the area; n – is the number of species of trees in the area.

The biomass stocks of tree trunks on each of the five sample plots were determined according to formula 4:

where: AGB_{g, j} – biomass stock of stem-wood calculated based on ground-data for the j-th sample plot, kg; n_j – number of trees on the j-th sample plot; ρ_m – average (across all plots) wood density, kg/m³; V_k,j – volume of the k-th tree on the j-th sample plot, m³.

The dependence of tree trunk diameters on height was determined using the multi-model technique. For this, the parameters of three simple regression allometric models were identified based on field data obtained from all five plots, and the median value of the results was calculated (equation 5).

$$d_{\mathrm{mod}}=\mathrm{median}(a\text{⋅}\exp (b\text{⋅}h),c\text{⋅}{h}^d,e\text{⋅}h),$$

(5)

where: d_mod – modeled trunk diameter (at a height of 1.3 m), m; h – tree height, m; model parameters: a, m; b, 1/m; c, m (this parameter includes the normalization factor h₀, m); e и d – dimensionless model parameters. The error in parameter identification was calculated using stochastic modeling [33,34]: random numbers with a uniform distribution in the range of five times the ground measurement error were added to the input data, which for trunk diameter ranged from ─0.05 to 0.05 m, and for tree height – from ─0.5 to 0.5 m); based on the resulting data, parameter identification of the model was performed by least-squares fitting (minimizing the sum of the squares of the differences between the measured and modeled diameter). This procedure was repeated 1000 times, resulting in 1000 sets of model parameters. For each parameter, the mean and standard deviation were calculated from its 1000 realizations. Models with identified parameters were then used to calculate the trunk volume based on their height obtained using lidar surveying, as described below.

2.2. LiDAR Data Acquisition and Processing

Lidar surveying was conducted in July 2023 in a mixed forest at the "Mukhrino" study site, covering an area of 662 hectares, using a UAV DJI Matrice 300 and a DJI Zenmuse L1 lidar. The UAV flight altitude was 120 meters above the surface, with the lidar oriented vertically downwards (nadir), and the scanning angle varied between ─35° and 35°. The point cloud density of laser reflections (LiDAR) was 102 points/m², ground sampling distance – 3.27 cm/pixel, flight speed – 12 m/s, transverse overlapping of LiDAR – 50%. Initial processing of the point cloud of laser reflections (LiDAR) included the alignment and smoothing of LiDAR to eliminate artifacts caused by distance measurement errors, noise, and UAV movements using DJI Terra, and then the data were exported in .las format.

Segmentation of individual tree crowns and calculation of their height based on LiDAR was performed using the "Lidar toolbox" module in MATLAB 2023a. The module requires the specification of two main parameters: "gridRes" – the size of the moving window for creating the canopy height model and "minTreeHeight" – the minimum tree height considered. In practice, for low-density LiDAR (about hundreds of points per m²), obtained for a forest with high crown closure and little variability in tree heights, the "Lidar toolbox" tends to underestimate the number of trees due to the "merging" of several plants into one during segmentation [19,31]. However, the degree of "underestimation" decreases as the "gridRes" parameter value decreases. Therefore, to analyze tree heights across the entire study site with the given LiDAR density, we selected the minimum functional value of the "gridRes" parameter = 0.3 m, and set the "minTreeHeight" to 5 m.

The result of the LiDAR segmentation was exported from MATLAB 2023a as two vector files in .shp format: a point file representing the tree crown vertex locations (with coordinates in the WGS84 42N projection in the attribute table) and a polygon file corresponding to the boundaries of segmented tree crowns.

2.3. Calculation of Carbon Storage And Their Ground Verification

As described above, the assessment of carbon stock in trunk wood was based on field data on the ratio of tree heights and trunk diameters (calculated using an ensemble of simple allometric models), the number of trees in the studied plots, and the wood density in them. The trunk diameters (d_mod) were calculated at the study site based on the height data obtained using lidar surveying (h_LIDAR), and from this, the cross-sectional area was calculated using equation 6:

$$s_{\mathrm{mod}}={\displaystyle \frac{\mathsf{\pi }\text{⋅}{d}_{\mathrm{mod}}^2}{4}},$$

(6)

where: s_mod – cross-sectional area of trees (at a height of 1.3 m), m²; d_mod – modeled trunk diameter (at a height of 1.3 m), m.

Based on the known height and cross-sectional area of trees, the merchantable volume (assuming the trunk has the shape of a cone) was calculated using equation 7:

$$V_m=h_{\mathrm{LIDAR}}\text{⋅}{s}_{\mathrm{mod}}\text{⋅}{{\displaystyle \frac{4}{3}}\text{⋅}\left({\displaystyle \frac{d_{1/2}}{d_{\mathrm{mod}}}}\right)}^2,$$

(7)

where: V_m – merchantable volume, m³; h_LIDAR – tree height obtained from lidar surveying, m; s_mod – cross-sectional area of trees (at a height of 1.3 m), m²; d_mod – modeled trunk diameter (at a height of 1.3 m), m; d_1/2 – diameter of the trunk at half the tree height (m), which was calculated using equation 8:

where: d_1/2 – diameter of the trunk (at half the tree height, m); h_LIDAR – tree height obtained from lidar surveying, m. The biomass stocks of tree trunks on each of the five sample plots were determined using equation 9:

where: AGB_{LIDAR, j} – biomass stocks of stem wood calculated based on LIDAR data for the j-th sample plot, kg; V_m,k,j – merchantable volume calculated using equation 7 for the k-th tree on the j-th sample plot, m³. The obtained merchantable volume of trunk wood (V_m) was converted to carbon stock using equation 10:

$$C=\mathsf{\alpha }\text{⋅}\mathrm{AG}{B}_{\mathrm{LIDAR}}\text{⋅}{C}_{\mathrm{frac}}$$

(10)

where: С – carbon stock in tree trunks, kg; α – dimensionless ground verification coefficient, identified by minimizing the sum of squares of the differences between the aboveground biomass stocks determined for each of the five sample plots both on the ground (AGB_g,i) and remotely (AGB_LIDAR,i); AGB_LIDAR – biomass stocks of tree trunks for the entire study site; C_frac – average weighted carbon fraction in tree trunks, calculated taking into account the number of coniferous and deciduous trees on the plots (similarly to ρ_w, according to the data [35] and [36]), and equal to 0.5.

3. Results

3.1. Forest Taxation and Allometric Models

In the result of ground forest inventory on sample plots 1–5, there were found 90, 66, 83, 92, and 59 trees in the overstory layer and 277, 196, 190, 202, and 212 in the understory layer respectively. The median (1Q, 3Q) height of trees in the overstory layer was 27 (24, 29), 26 (17, 28), 29 (23, 31), 27 (25, 29), and 29 (27, 30) m, and in the understory layer it was 15 (14, 17), 14 (13, 15), 14 (11, 17), 15 (13, 18), and 15 (13, 18) m. The median (1Q, 3Q) diameters of tree trunks in the overstory layer at a height of 1.3 m were 29 (24, 34), 25 (16, 33), 34 (24, 42), 29 (24, 36), and 41 (29, 49) cm, and in the understory layer they were 13 (11, 14), 12 (11, 13), 14 (10, 16), 13 (12, 16), and 13 (12, 16) cm. The average weighted wood density ρ_w varied slightly across the sample plots and was on average 0.43±0.02 t/m³. The biomass stocks AGB_g,i were 152, 107, 156, 148, and 137 t/ha respectively.

The parameters of allometric models (equation 5) for the dependence of tree diameter on height, identified using ground data (Table 2), as well as the type of parameterized model (Figure 2) are presented below.

3.2. Characteristics of the Forest stand Based on LiDAR data

In the result of assessing the stand characteristics based on LiDAR data, 41, 30, 37, 44, and 27 trees were found on sample plots 1–5 respectively. The number of trees identified using the LiDAR survey was significantly underestimated compared to the results of ground forest inventory (Figure 3): both for the number of overstory trees and for overstory and understory trees combined. Nevertheless, an extremely strong linear regression relationship (R² = 0.98) was found between the number of overstory trees determined on the ground and based on the LiDAR survey, while the relationship between the total number of trees determined on the ground and remotely was relatively low (R² = 0.36).

The median (1Q, 3Q) height of trees (without layer separation) obtained based on LiDAR data was 28 (26, 29), 27 (25, 29), 30 (29, 31), 28 (25, 29), and 28 (26, 30) m. The Mann-Whitney U-test, conducted to compare the heights obtained from LiDAR data and ground assessment (with separation into all trees and overstory trees) showed no significant differences only between LiDAR and overstory trees on plots 4 and 5 (p = 0.29 and 0.75 respectively, Figure 4). It is evident that LiDAR data with the survey parameters specified in section 2.2 are insufficient for inventorying trees below 22–24 m in height, which does not allow for height assessment beyond the overstory layer. In most cases (sites 1–4), the median height of trees obtained from LiDAR data is 1–2 m higher than that obtained using ground data for the overstory layer.

The modeled median (1Q, 3Q) diameter of tree trunks (without layer separation) at a height of 1.3 m was 33 (29, 39), 32 (26, 38), 39 (34, 44), 33 (27, 38), and 34 (29, 41) cm (Figure 5). The uncertainty (1Q, 3Q) in the results of modeling the median tree diameter on the plots includes two types of errors caused by: 1) natural variability in tree diameters on each plot and 2) error in model parameter identification; 1Q was calculated for the results of d_mod obtained from an ensemble of models using parameters identified at the lower uncertainty bound, and 3Q was calculated at the upper uncertainty bound (Table 2). As with the assessment of tree heights, the medians of the modeled diameters based on LiDAR data significantly differed from the ground assessment (except for site 3, p = 0.28). In most cases (sites 1–4), the median diameter of trees obtained from LiDAR data was 4–7 cm larger than those obtained using ground data for the overstory layer.

The modeled total biomass stocks for plots 1–5 (AGB_LIDAR,i) thus included the combined uncertainty of modeling the median diameter and assessing the tree height, and amounted to 65 (55, 75), 44 (38, 50), 80 (65, 95), 71 (59, 83), and 46 (38, 54) t/ha respectively. Nevertheless, they closely correlated with AGB_g,i obtained from ground data (Figure 6).

3.3. Ground verification of AGB and Calculation of Carbon Storage

The underestimation of biomass stocks obtained from LiDAR data was corrected by introducing a correction factor α = 2.2 (equation 10), used to minimize the RMSE. The results of the biomass stock adjustment are presented in Figure 7: RMSE was reduced from 80 to 32, approaching a 1:1 ratio between ground and remote assessments.

The obtained estimate αAGB_LIDAR was used to calculate the carbon stocks for the entire study area (S = 662 ha), which amounted to 35814±6915 tC or 54±11 t/ha. The spatial variability of carbon stocks in the study area is presented in Figure 8: it ranged from 5±1 tC/ha in areas bordering wetlands to 110±21 tC/ha in forest areas occupied by cedar and mature mixed stands.

4. Discussion

The calculated biomass stocks before the introduction of the correction factor α differed by 2.2 times from those obtained through ground assessment. We assume that the main contribution to this underestimation came from the unrecorded overstory and understory trees during the LiDAR survey: the number of identified overstory layer was underestimated by 2.2 times, and the understory layer of trees were completely ignored (as confirmed by the median tree height assessment on the trial plots – Figure 4). The reason for ignoring the understory trees is the high canopy closure of the overstory layer, which were themselves undercounted due to the low density of the point cloud during the LiDAR survey. To obtain a higher point cloud density, it is necessary to increase lateral overlap during the LiDAR survey, reduce the flight altitude, or conduct multiple surveys of the same area. Accounting for lower layer trees could be done in the fall period with fallen leaves, which could increase the number of reflected LiDAR signals from beneath the forest canopy. Nevertheless, the high correlation of tree counts obtained by ground (for the overstory layer) and remote assessment indicates a systematic error that can be effectively compensated, as shown above.

The systematic overestimation of the median tree height is insignificant (1–2 m, except for plot 5) and can also be explained by the unaccounted trees below 22–24 meters, which might have shifted the median assessment. The overestimation of tree diameters calculated (by 4–7 cm, again except for plot 5), may be due to the assumption that trees take the shape of a cone. It is noteworthy that the absence of diameter and height overestimation is characteristic only for plot 5, where the least number of trees were found among all the plots, suggesting that the process of obtaining reliable LiDAR surveys and tree identification is more accurate for sparse stands.

We compared the carbon stocks of the forest we studied with estimates available in the literature (Table 3). In the table, we showed the results of the current study for forest plots on relatively drained sites, and it largely coincides with the literature data for similar forests. The average carbon stock of stem wood across the entire study area is close to the lower limit of estimates for similar forests, due to the presence of a large number of bog areas (about 9% of the territory) as well as a river crossing the forest and its floodplain (7%). Of the 662 hectares, the forest areas directly occupy about 556 hectares, many of which are represented by low-productivity communities: if the total carbon stock is recalculated for the actual area they occupy, the specific carbon stock will be 64 t/ha. At the same time, the average value obtained from all analyzed works is 83 t/ha. It should be noted that in the latter case some authors included not only stem wood but also the biomass of branches and photosynthesizing organs of trees, which was not the aim of our study.

Testing the technology for accounting the total aboveground biomass of trees, as well as accounting for the spatial variability of their species composition, will be the topic of our further work.