Mapping and Estimating Aboveground Biomass in an Alpine Treeline Ecotone under Model-Based Inference

Abstract

Due to climate change, treelines are moving to higher elevations and latitudes. The estimation of biomass of trees and shrubs advancing into alpine areas is necessary for carbon reporting. Remotely sensed (RS) data have previously been utilised extensively for the estimation of forest variables such as tree height, volume, basal area, and aboveground biomass (AGB) in various forest types. Model-based inference is found to be efficient for the estimation of forest attributes using auxiliary RS data, and this study focused on testing model-based estimations of AGB in the treeline ecotone using an area-based approach. Shrubs (Salix spp., Betula nana) and trees (Betula pubescens ssp. czerepanovii, Sorbus aucuparia, Populus tremula, Pinus sylvestris, Picea abies) with heights up to about five meters constituted the AGB components. The study was carried out in a treeline ecotone in Hol, southern Norway, using field plots and point cloud data obtained from airborne laser scanning (ALS) and digital aerial photogrammetry (DAP). The field data were acquired for two different strata: tall and short vegetation. Two separate models for predicting the AGB were constructed for each stratum based on metrics calculated from ALS and DAP point clouds, respectively. From the stratified predictions, mean AGB was estimated for the entire study area. Despite the prediction models showing a weak fit, as indicated by their R²-values, the 95% CIs were relatively narrow, indicating adequate precision of the AGB estimates. No significant difference was found between the mean AGB estimates for the ALS and DAP models for either of the strata. Our results imply that RS data from ALS and DAP can be used for the estimation of AGB in treeline ecotones.

1. Introduction

Forests play a major role globally as carbon sinks. Hence, afforestation may be important for climate change mitigation through carbon sequestration [1]. However, there are large differences in the capacity of different forest types to sequester carbon and produce biomass. These differences depend on key environmental factors such as the available nutrients in the soils and the climatic conditions. In the Nordic countries, the lowland and low-latitude forests constitute the most productive areas and the largest pools of tree biomass, but even high-elevation ecosystems are known to store large amounts of carbon in vegetation and the soil [2,3,4].

Climate changes affect the establishment, growth, and mortality of trees and other woody vegetation [5,6,7,8], and the effects might be more pronounced in the transition zones between the boreal and alpine zone (i.e., the treeline ecotone), where trees grow at their tolerance limit in terms of temperature compared to productive forests at lower elevations [9,10]. Therefore, in the treeline ecotone, even a moderate increase in temperature or change in precipitation might lead to the increased growth of existing trees and promote the establishment of pioneer trees in currently treeless areas [11,12]. According to Bryn et al. [13], the alpine areas in Norway where pioneer trees potentially could establish constitutes between 25 and 30% of the total land area (excluding bare rock and barren areas). These are potential areas for forest expansion where the potential impact of carbon sequestration is substantial. An increase in tree vegetation into the current alpine areas will also have an impact on albedo [14,15,16], especially during winter as previously entirely white surfaces will have scattered dark areas which absorb more incoming solar radiation.

There are several other factors apart from climate change that could induce changes in the treeline ecotone, such as herbivory [17,18]. In Norway, where the current study was carried out, it has been common in the past for domestic animals to graze and browse in montane areas during the summer months (summer farming), in some areas more intensively than others. Considering the large areas in play, the decline in summer farming and herbivory can have a substantial effect on the biomass and carbon stocks in high-elevation forests. Previous research seeking to disentangle and quantify the importance of the different factors affecting the treeline ecotone has highlighted the complexity of the causal relationships of treeline dynamics [5,17,18,19]. Even though both climate and herbivory have been suggested as important factors, they only explain a small proportion of the variation in the observed changes in the treeline ecotone. This likely means that the climatic and herbivory variables fail to fully represent the true effects of these factors and that there are context-dependent local factors that go unaccounted for in such analyses [9,10,20,21,22,23]. Since climate-induced responses in the treeline ecotone are difficult to predict, the future development of biomass stock is highly uncertain. A monitoring system where objective data are collected on a regular basis that enables the estimation of actual biomass and carbon stocks is therefore important to fulfil national obligations with regard to carbon reporting.

Accurately estimating and monitoring changes in forest biomass and carbon content is crucial for meeting the requirements of both the Kyoto Protocol and the Paris Agreement [24]. In Norway, the national forest inventory (NFI) is the primary source for biomass and carbon estimates. For forest areas, the sampling grid of the NFI is 3 km × 3 km, but for montane regions, it is 3 km × 9 km [25], which is sparse in terms of providing precise local and overall country-wide estimates. However, with the use of auxiliary remotely sensed (RS) data, the acquisition of data with wall-to-wall coverage for large regions has, in many cases, been shown to be effective in terms of obtaining precise estimates of variables of interest [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41]. Compared to pure field-based samples, RS data offer the advantage of quickly providing coverage over large areas, and wall-to-wall RS data are often available over the entire area of interest (AOI). Another advantage of RS data is their ease of use in remote or inaccessible areas, which is often the case where treeline ecotones are found.

Numerous studies have demonstrated the versatility of three-dimensional (3D) point cloud data obtained through airborne laser scanning (ALS) for estimating the properties of forest and vegetation (e.g., [34,37,39,42,43,44,45,46,47]). Næsset and Nelson [48] found that almost all treeline trees with a height of 1 m or taller could be detected using ALS data. However, trees below 1 m in height were often not detected because a laser pulse needs to intercept a minimum surface area before an echo is triggered. With point data, the probability of detecting trees will also depend on the point density. They also found that ground features like rocks could yield echoes with positive height values because they were misclassified as vegetation echoes. Thus, although detecting the smallest trees in the treeline ecotone using ALS presents challenges, there is potential for developing an effective monitoring system for AGB based on ALS, particularly for trees and shrubs taller than 1 m. Despite the steep decline in detection rates of single trees shorter than 1 m, an area-based approach may still be a viable option for estimating AGB. Some studies have already used bi- and multi-temporal ALS for estimating changes in height among treeline trees [49,50], and the results from these studies are encouraging, showing the potential utility of ALS for monitoring purposes.

The use of 3D point cloud data from digital aerial photogrammetry (DAP) is an alternative to ALS, especially because of the potential finer spatial resolution. Aerial images are cheaper to acquire and can therefore be more frequently collected compared to ALS data [31]. Therefore, DAP data has emerged as an alternative to ALS in operational forest inventories when both cost and utility are taken into account [30,51,52]. Studies that have applied DAP for detection and estimation purposes related to small trees are still few, although Puliti et al. [53] studied the use of aerial images from a drone platform for the estimation of biophysical properties in productive forests under regeneration. In the latter study, it was found that the height values of DAP point clouds tended to underestimate tree height more than those from ALS. It was also found that solitary single trees were more likely to be smoothed from the DAP point cloud. Whether or not these results were directly related to the settings and algorithms used in the matching of the images was not identified. Therefore, further investigation into the use of DAP for vegetation attribute estimations in such transitional zones was recommended.

The uncertainties associated with estimates of state and changes in the AGB using auxiliary RS data have been studied for mature forest stands [29,54,55,56,57,58,59,60], montane forests [46], and young forests [61]. The estimation of AGB using auxiliary RS data can be accomplished through design-based or model-based inferences, which have been discussed in detail in Ståhl et al. [33]. The model-based inference is independent of the probability sampling assumptions, which makes it advantageous for inaccessible and remote study sites. However, the accuracy of the estimates is solely dependent on the applied model. Thus, correctly specified models are important to avoid systematic model prediction errors [33,43,62,63]. For informed decision making, it is also important that estimates of precision are provided along with the estimates of AGB. In a study by Næsset et al. [31], the height of treeline vegetation was estimated, and the efficiency of ALS and DAP RS data was compared by characterising the uncertainty in the height estimations using root mean squared error. However, to the best of our knowledge, there have been no studies where the uncertainty of AGB estimates in treeline ecotone sites has been estimated. Therefore, further research is needed to better understand the uncertainties associated with estimating AGB in this region.

Objective

In the current study, AGB for a treeline ecotone site was estimated with corresponding estimates of precision using an area-based approach and model-based inference. Models for field measured AGB were constructed using metrics calculated from both ALS and DAP data. The study area was divided into two strata (tall and short vegetation), and separate models were fitted for each stratum with both ALS and DAP metrics. Along with the estimates of total AGB based on both sets of metrics, standard error estimates were obtained by means of parametric bootstrapping (PB). The main objective was to evaluate and compare the precision of AGB estimates obtained utilising ALS and DAP as auxiliary data and to assess the possibilities of using the area-based approach aided by these RS data sources to estimate AGB in treeline ecotones.

2. Materials and Methods

2.1. Study Area

The AOI (7°58′E, 60°42′N; Figure 1) is located in the municipality of Hol, Norway, at elevations between 1050 and 1320 m above sea level [21]. The dominant treeline species in this region is mountain birch (Betula pubescens ssp. czerepanovii), along with a few individuals of rowan (Sorbus aucuparia), aspen (Populus tremula), Scots pine (Pinus sylvestris), and Norway spruce (Picea abies). The shrub species considered in this study were Salix lapponum, Salix glauca subsp. glauca, and Betula nana. Other woody dwarf shrubs, such as Vaccinium spp. and Empetrum sp., are present in the area but not considered in the current study. For an overview of all plant species and their abundances within the study area, see [64].

A graphical overview of the workflow is displayed in the chart in Figure 2. The data collection, modelling, and estimation procedures are discussed further in the following sections.

2.2. Field Methods

There are both destructive and non-destructive methods available for measuring the aboveground biomass (AGB) of trees; they differ in terms of resource requirements and operational considerations [65,66]. In this study, our decision to employ a non-destructive method for measuring AGB values in the field was primarily driven by resource limitations. Destructive sampling, while providing more precise biomass estimates, would have required significant additional resources that were not available to us. Furthermore, in an operational setting, where resource allocation is constrained, the use of destructive sampling would not be feasible or practical. Lastly, it is worth noting that, in this specific case, we also faced restrictions and did not have permission to perform destructive sampling. As a result, we utilised a predictive approach, leveraging field measurements to estimate the total AGB for each plot. Despite the inherent limitations, this approach enabled us to derive biomass estimates effectively while avoiding the challenges associated with destructive sampling. The field data were collected partly using wall-to-wall ALS data to guide the location of the sample plots. The area was tessellated into 100 m² grid cells, and the average height ($H_{mean}$) of first echoes of the ALS pulses within each grid cell was calculated. Based on the $H_{mean}$ values, the study area was split into two strata. Grid cells where $H_{mean}$ < 1 m were categorised as ‘short woody vegetation’ (SHORT), and grid cells where $H_{mean}$ ≥ 1 m were categorised as ‘tall vegetation’ (TALL). After classification into SHORT and TALL, the SHORT grid cells were further tessellated into 16 equally sized squares (6.25 m²).

Twenty TALL grid cells with $H_{mean}$ values evenly distributed over the entire range of $H_{mean}$ values were selected to establish field reference plots to initiate modelling. The selection process involved dividing the $H_{mean}$ range between 1 m and the maximum value into 10 equally sized bins. Then, we randomly selected two grid cells from each bin. The fieldwork was carried out in summer 2019. A Topcon HiPer SR geodetic-grade GNSS receiver in real-time kinematic (RTK) mode was used to navigate to the centre of each selected grid cell, where a circular plot with a radius of 5.64 m (area = 100 m²) was established. Within each plot, all trees with diameter at breast height (dbh) > 0 were calipered. Breast height was defined at 1.3 m above ground. Tree height was measured using a folding rule or a Vertex hypsometer on two subjectively selected height-sample trees in each plot. The selection was made so that the 40 selected trees (two trees per plot) covered a range of dbh values. In addition, a 1.5 m radius (7.07 m²) sub-plot was established at the centre of each plot to sample shrubs and tree saplings with a height (h) < 1.3 m. All tree species were sampled, including juniper (Juniperus communis), as well as the shrubs species willow (Salix spp.) and dwarf birch (Betula nana). The sub-plot was divided into four quadrants by two perpendicular lines that intersected at the plot centre in north–south and east–west directions (Figure 3). In each quadrant, the tree sapling or shrub closest to each of the points 1 m from the plot centre in the directions that were 45, 135, 225, and 315 degrees relative to north (Figure 3) was measured for height and diameter at root collar. In the southeast and northwest quadrants, the number of saplings and shrubs were counted.

The field measurement for the SHORT stratum took advantage of 180 previously established vegetation monitoring plots [20]. The original vegetation plots were squares of 0.25 m², but we established 1.5 m radius plots in the NW corner of each square. Two additional plots were purposely established in tall Salix vegetation, since this was poorly covered by the 180 prepositioned plots. The same measurements as those described for the sub-plots of the 20 plots of the tall vegetation were carried out. The position at the centre of each plot was registered with the same Topcon receiver described above in RTK-mode.

For each plot, total AGB was predicted as the sum of individual tree AGB, following the same procedure as in [67,68]. For the TALL plots, the predictions were separate for trees with h < 1.3 m and those with h $\ge$ 1.3 m. For trees with h $\ge$ 1.3 m, the base heights were first predicted using the base height models of Fitje et al. [69]. Then, single-tree biomass models [70] with field-measured dbh and predicted base height as inputs were used to predict the base biomass values (b1) of each tree. For each of the height-sample trees, a second biomass (b2) value was predicted using measured dbh and measured h. Then, a common correction factor k for the base-biomass was calculated from all the height sample trees as the ratio between the sum of b2 and the sum of b1. Biomass for each tree was predicted as b1 × k. Single-tree biomass predictions were summed for each plot and scaled to a per hectare value. For trees and shrubs with h < 1.3 m observed on the 1.5 m radius sub-plots, a mean diameter at the root collar and a mean height were calculated from the measurements in each plot. These mean values were used as inputs in the single-tree biomass model of Kolstad et al. [71] to predict mean-tree biomass and then multiplied with the number of individuals and scaled to a per hectare value. The predicted biomass for trees with h < 1.3 m was added to the biomass prediction for the taller trees. For the SHORT vegetation plots, the total biomass was calculated using the same procedure used for the 1.5 m radius sub-plots of the TALL vegetation stratum. A summary of field reference data is presented in

Table 1. Summary of field reference data collected in 2019.

Stratum	n	Mean AGB (Mg ha−1)	SD (Mg ha−1)	Min (Mg ha−1)	Max (Mg ha−1)
TALL	20	30.2	17.5	8.48	65.8
SHORT	182	1.61	2.68	0.00	18.3

.

2.3. Remotely Sensed Data

2.3.1. Data Acquisition and Initial Processing

The sensor and flight information for the RS data is provided in

Table 2. Summary of sensor and flight information for ALS and DAP.

	ALS	DAP
Sensor system	Riegl VQ-1560i	Sensefly S.O.D.A. camera
Platform	Piper PA-31-350 Chieftain	Sensefly eBee
Acquisition dates	8 and 25 June 2018	7–10 July 2019
Flight altitude (m a.g.l) *	3400	120
Flight speed (m s−1)	NA	15
Point repetition frequency (KHz)	350	NA
Scan frequency (Hz)	162	NA
Point density (points m−2)	2	55
Half scan angle (degrees)	20	NA

* above ground level.

. ALS data were collected by the contractor, Terratec AS, as part of national scanning, with a point density of 2 points/m². The raw point clouds were pre-processed by the contractor, and laser echoes were classified as “ground”, “unclassified”, “noise”, “bridge”, or “snow”. Planimetric coordinates and orthometric heights were computed for all echoes. A triangulated irregular network (TIN) was then generated from the laser echoes classified as “ground” using the R package “lidR” [72,73]. Since the study aimed to provide estimates of the biomass of shrubs and trees whose lower stems often grow parallel to the ground due to snow pressure, all vegetation echoes were included without setting a threshold on echo height. Thresholds of <2 m have commonly been used in previous studies to omit falsely classified vegetation echoes [50,74,75].

The DAP data were collected following a flight plan with perpendicular flight lines with lateral and longitudinal overlaps of 80%. Across the entire study area, 43 orange wooden crosses (30 cm width) were distributed uniformly and positioned and used as ground control points (GCP) to enable computation of the position and orientation of the images. The DAP point cloud was thinned from a mean density of 55 points m⁻² to 38 points m⁻² using the R package “lidR” [72,73] to obtain a uniform point density over the entire area. Examples of ALS and DAP point clouds are presented in Figure 4 for two grid cells, one in each of the TALL and the SHORT strata. The higher point density of the DAP point clouds is clearly visible.

2.3.2. Correction of the DAP Point Cloud

An initial inspection of the normalised DAP point cloud revealed many points with negative normalised height values. Regions with no observable vegetation, meaning they should have heights of 0 m, were observed to contain points with both positive and negative heights. These were observed to be in the magnitude of 0–2 m. This indicated that the error was not simply a global vertical shift and that this error was too large to be solely due to the inaccuracy of the DAP system. The most likely culprit was the positioning accuracy of the GCP markers. Positioning errors could result in either rotation or translation in the DAP registration and thus errors in both the horizontal and vertical directions.

The ALS and DAP point clouds did not share any GCP markers. Furthermore, the density of the ALS data made the extraction of natural terrain features that could be applied as GCPs too inaccurate. As a result, the DAP point cloud registration could not be corrected by simply using shared tie points in the ALS and DAP point clouds. As an alternative to this approach, the local vertical errors between the point clouds were estimated and used to correct the original DAP point cloud. This approach can be summarised in three steps: (i) ground points in the DAP point cloud were classified using the software TerraScan [76], (ii) the ALS ground point heights were interpolated at the XY coordinates of the DAP ground points so that the vertical shift between ALS and DAP ground heights could be computed, and (iii) the correction of the entire DAP point cloud was interpolated from the results of the previous step. The interpolation was performed using the “scipy.interpolation.griddata” Python package. If the point to be interpolated was located inside the convex hull of the input points, linear interpolation was used; otherwise, nearest-neighbour interpolation was employed.

2.3.3. Computation of Metrics

For both the ALS and DAP point clouds, metrics that represented the height distribution and the density of the point cloud for each plot and grid cell were computed using the echoes classified as “first of many” and “single”. The normalised ALS and DAP point clouds were clipped to the spatial extent of the plots and grid cells, and the points with negative heights were removed from further processing for the TALL stratum. For the SHORT stratum, the DAP points with negative heights were kept to account for small height shifts still present in the point cloud after correction. We considered this to be important because, for the SHORT stratum, most of the AGB was growing close to the ground. Removing negative heights could possibly equalise the variation in the DAP point cloud in areas where AGB actually varies. The metrics representing height were heights at the deciles of the echo height distribution. The densities were the ratios between the number of echoes above 10 different height-levels and the total number of echoes. In addition, the standard deviation of the heights, skewness, and kurtosis were computed; all together, 23 metrics for each of the ALS and DAP point clouds were considered.

2.4. Model Construction

Since AGB was small within our study area, the values of both the ALS and DAP metrics were zero for many of the field plots, which was a limitation concerning variable transformations and model forms. After some preliminary tests with different modelling approaches, such as zero-inflated and log-transformation of the response variable, we chose to construct linear models fitted using ordinary least squares. Linear models were suitable in this case because the training dataset contained extreme AGB values, which reduced the need for extrapolation. The models were constructed separately for the TALL and SHORT strata. For the TALL stratum, the explanatory variables for the ALS and DAP datasets were separately selected through a stepwise procedure based on the Bayesian information criterion (BIC) and adj-R² [77]. The general form of the models for the TALL stratum was as follows:

$$\mathrm{AGB}={\beta }_0+{\beta }_1{X}_1+\cdots +{\beta }_j{X}_j+\epsilon$$

(1)

where j is the number of explanatory variables (X), ${\beta }_0$, ${\beta }_1$, …, ${\beta }_j$ are the parameter estimates, and $\epsilon$ is the random error term.

As indicated in Section 2.2, there was a difference in the area of the field plots for the short vegetation (7.07 m²) and the area of the grid cells (6.25 m²). Thus, a scale-independent explanatory variable was purposely chosen for the models for the SHORT stratum. The metric chosen for both the ALS and DAP models was the mean point height ($H_{mean}$), which was retained even if the parameter estimate was statistically non-significant.

For the SHORT stratum, the models were formulated as follows:

$$\mathrm{AGB}={\beta }_0+{\beta }_1{H}_{mean}+\epsilon$$

(2)

where $H_{mean}$ is the average height of either the ALS or DAP points, β₀ and β₁ are the model coefficients, and $\epsilon$ is the random error term.

The models’ training accuracies were assessed using the root mean square residual error (RMSE) and relative RMSE (rel.RMSE) [78]:

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\widehat{\mathrm{AGB}}}_i-{\mathrm{AGB}}_i\right)}^2}{n}}$$

(3)

$$\mathrm{rel}.\mathrm{RMSE}=\frac{\mathrm{RMSE}}{{{\displaystyle \sum }}_{i=1}^n\left({\mathrm{AGB}}_i\right)/n}\times 100$$

(4)

where ${\widehat{\mathrm{AGB}}}_i$ is the model-predicted AGB for field plot i, ${\mathrm{AGB}}_i$ is the corresponding field reference AGB, and n is the number of observations in the reference data used to train the models.

2.5. Estimation of Mean AGB

The final AGB—ALS and AGB—DAP models were applied to the grid cells belonging to the stratum for which they were constructed. Mean AGB was estimated separately for both strata and with both models (${\overline{\widehat{\mathrm{AGB}}}}_{TALL}$, ${\overline{\widehat{\mathrm{AGB}}}}_{SHORT}$) as the means of the respective stratified AGB predictions. Mean AGB values for the entire area (${\overline{\widehat{\mathrm{AGB}}}}_H$) using both models were estimated as area weighted means for both strata individually for ALS and DAP as [79]:

$${\overline{\widehat{\mathrm{AGB}}}}_H=\frac{{\overline{\widehat{\mathrm{AGB}}}}_{TALL}. A_{TALL}+{\overline{\widehat{\mathrm{AGB}}}}_{SHORT}. A_{SHORT}}{A_H}$$

(5)

where A is the area, and the subscripts TALL, SHORT, and H denote the TALL stratum, SHORT stratum, and the entire area, respectively.

2.6. Variance Estimation via Parametric Bootstrapping

Estimation of the standard error (se) for the AGB estimates was carried out using parametric bootstrapping (PB). PB is based on Monte Carlo simulation and is convenient within the model-based inference framework [80]. It enables the statistical inference of the variable of interest when the true distribution is unknown [81], and the Monte Carlo errors become negligible when the PB samples are large [46]. In this study, the number of bootstrap samples (${\mathrm{N}}_{\mathrm{PB}}$) was set to 50,000, and for each iteration (k), the mean of AGB predictions (${\overline{\widehat{\mathrm{AGB}}}}_k$) was calculated. The standard error for the AGB predictions ($\widehat{se}$) was obtained separately for each stratum as [46]:

$$\widehat{se}=\sqrt{\frac{1}{{\mathrm{N}}_{\mathrm{PB}}-1} {\displaystyle \sum }_{k=1}^{{\mathrm{N}}_{\mathrm{PB}}}{\left({\overline{\widehat{\mathrm{AGB}}}}_k-{\overline{\widehat{\mathrm{AGB}}}}_{\mathrm{PB}}\right)}^2}$$

(6)

where ${\overline{\widehat{\mathrm{AGB}}}}_{\mathrm{PB}}$ is the mean of the ${\overline{\widehat{\mathrm{AGB}}}}_k$ values of the ${\mathrm{N}}_{\mathrm{PB}}$ samples. To assess whether the number of iterations was sufficient to accurately estimate the uncertainty of AGB predictions, we calculated a standard error stabilisation indicator (se.s). This indicator evaluates the change in the standard deviation of mean predicted AGB as additional bootstrap iterations are conducted. The se.s value was determined after each iteration as the maximum difference between the standard deviation of mean predicted AGBs obtained during the last 50% of iterations and the standard deviation of all mean predicted AGBs obtained in the final iteration. According to the authors of [80], an se.s of <0.5% of the se after the last iteration indicates that the se has stabilised sufficiently.

The area-weighted standard error for the entire study area $\left({\widehat{se}}_H\right)$ combining both stratum estimates can be estimated as [79]:

$$\widehat{{\mathrm{se}}_H}=\sqrt{\frac{\left({\widehat{\mathrm{se}}}_{TALL}^2.A_{TALL}^2+{\widehat{\mathrm{se}}}_{SHORT}^2.A_{SHORT}^2\right)}{A_H^2}}$$

(7)

The 95% confidence intervals (CI) were obtained for the mean estimate of predicted AGB for the respective strata and the total area as [82]:

$$\mathrm{CI}=\overline{\widehat{\mathrm{AGB}}} \pm \left(\widehat{se}\times 1.96\right)$$

(8)

$${\mathrm{CI}}_H={\overline{\widehat{\mathrm{AGB}}}}_H \pm \left({\widehat{se}}_H\times 1.96\right)$$

(9)

3. Results

For both of the TALL stratum models, the 9th height decile ($H_{90}$) was selected as the only explanatory variable, while for the SHORT stratum, the scale-independent $H_{mean}$ was selected. The estimated parameter coefficients for all models are displayed in

Table 3. The table displays the explanatory variables, model coefficients, adj-R², RMSE, and relative RMSE (rel.RMSE) of the AGB prediction models for both strata and data sources.

Stratum	Model	Explanatory Variable	Prediction Model	adj-R2	RMSE (Mg ha−1)	rel.RMSE (%)
TALL	AGB-ALS	$H_{90}$	$-15.09+11.64H_{90}$	0.47	12.4	41.1
	AGB-DAP	$H_{90}$	$-0.388+8.25H_{90}$	0.43	12.8	43.8
SHORT	AGB-ALS	$H_{mean}$	$0.74+32.3H_{mean}$	0.15	2.47	154.2
	AGB-DAP	$H_{mean}$	$0.96+28.41H_{mean}$	0.27	2.28	118.1

along with the corresponding adjusted R² (adj-R²) and root mean squared error (RMSE). For all models, the adj-R² values were relatively small, indicating moderate to weak model fit. The relatively large proportion of unexplained AGB variation can also be seen in Figure 5, where field-observed AGB is plotted against the model fitted values for all four models. There were no large differences in model fit between the models based on ALS and DAP data, but substantial model fit differences between the strata were observed. Table 3 also shows that the slope parameter values for the prediction models are smaller for the DAP model for both strata compared to those of the ALS models, indicating that the two point clouds have different height distributions.

Figure 6 displays the se.s after each bootstrap iteration for both model types and both strata. In all data types and stratum combinations, the se.s dropped below 0.5% before reaching 50,000 iterations. However, in some instances and to a varying degree, se.s increased in value after the stabilisation criterion was first reached; however, at the final iteration of all four simulations, the se.s remained below the threshold of 0.5% for at least 10,000 of the last iterations.

Table 4. Estimation results of the AGB prediction models. $\overline{\widehat{\mathrm{AGB}}}$ is the estimated mean, $\widehat{\mathrm{se}}$ the corresponding standard error estimate, and $\widehat{\mathrm{CI}}$ the 95% confidence interval.

Stratum *	Model	$\overline{\widehat{\mathbf{A}\mathbf{G}\mathbf{B}}}$ (Mg ha−1)	$\widehat{\mathbf{s}\mathbf{e}}$ (Mg ha−1)	$\widehat{\mathbf{C}\mathbf{I}}$ (Mg ha−1)
TALL A = 0.97 ha	AGB-ALS	26.5	3.16	23.4–29.7
	AGB-DAP	29.2	3.13	22.9–35.5
SHORT A = 270.75 ha	AGB-ALS	2.05	0.20	1.66–2.45
	AGB-DAP	1.93	0.17	1.59–2.27
Total A = 271.72 ha	AGB-ALS	2.14	0.39	1.75–2.53
	AGB-DAP	2.03	0.34	1.68–2.37

* A = area in ha.

displays the results from the estimations for each stratum and model, as well as for the entire study area. There were no significant differences between the values of $\overline{\widehat{\mathrm{AGB}}}$ between models or within strata. Due to the large difference in the size of the area regarding the two strata, the area-weighted values of $\overline{\widehat{\mathrm{AGB}}}$ for the entire study area were similar to those for the SHORT stratum.

4. Discussion

This study employed a model-based approach with data from ALS and DAP as auxiliary information to estimate AGB and its corresponding precision. In contrast to other comparable studies [46,74], our AOI was exclusively covered by shrubs and short trees. As a result, this research has expanded the boundaries of RS technology in terms of supporting AGB estimation. Despite the small AGB values in our AOI, the widths of the CIs were encouraging in terms of operationalising AGB estimation in the treeline ecotone using both ALS and DAP as auxiliary information. The 95% CIs around the estimated mean AGBs overlapped for ALS and DAP in both strata (Table 4 and Figure 5), indicating that the estimates derived from ALS and DAP were not statistically significantly different. However, the relative uncertainties were larger compared to those obtained in previous studies for productive (e.g., [40,41,54,83]) and montane forests [46] due to the irregular shapes of stems and tree crowns in our AOI, which created weaker relationships between AGB and the RS data and led to higher relative measurement errors. Compared to productive and montane forests, where trees are substantially taller than in our AOI, the relative magnitude of the measurement errors of the RS data was also expected to be larger.

Preprocessing of the DAP point cloud was carried out to correct varying degrees of shifts in point heights or horizontal errors over the study area. Since our AOI comprised mostly short vegetation, this preprocessing was deemed highly important. The corrected values of the original DAP ground points were centred around 0, and both negative and positive corrections were carried out, indicating that the approach was at least partly successful at correcting erroneous height values at ground level. Errors still present after correction were considered to be random. A weakness of our approach was that it assumed vertical errors only. However, errors in the XY directions also must be expected, although they were most likely small. The impact of such potential horizontal errors would have been larger with a reduced size of the field reference plots used for modelling [84,85] and would have been manifested by larger uncertainties in the model parameters. The uncertainty will increase with increasing spatial heterogeneity of the vegetation over the AOI. In the current study, there is reason to believe that such errors had the greatest impact on the SHORT stratum since the field plot size was small. In the prediction phase, the horizontal errors will be less important since there is no direct requirement of spatial coherence. However, the precision of the predictions would still be affected by the errors associated with the model.

Judging by adj-R², the model relationships were weaker in the current study compared to studies carried out on boreal productive forests (e.g., [34]) mainly composed of conifer trees with regular conical crown shapes. Due to the irregularity and large variability of stem and crown shapes in our AOI, in addition to our population comprising mainly of shrub species and small individual trees, this difference was not surprising. In the SHORT stratum, most of the trees and shrubs were short and extended their stems partly along the ground, whereas the trees that constituted the TALL stratum were often snow-bent with stems that were both crooked and growing partly along the ground. As a result, a large part of the AGB was found within the uncertainty range of the TIN surface, which is often around 20–30 cm [86].

The positional errors of the plots would have had similar effects as those from the horizontal shift in the point cloud, as discussed above. However, in the open montane landscape of our AOI, the conditions for GNSS positioning were generally excellent, reducing the potential for positional errors. Nonetheless, potential positional errors are of greater importance for the SHORT stratum than the TALL stratum due to the smaller plot size used.

Furthermore, vegetation near the plot boundary can cause edge effects [85] as trees and shrubs rooted just outside the plots may intersect with the plots. In such cases, the trees or shrubs would not be measured in the field; therefore, they would not be part of the measured AGB, but they would still affect the properties of the point clouds for these plots. Conversely, trees and shrubs inside the plots close to the plot edge could lean out of the plot and still be included in the AGB estimate, affecting the properties of the point clouds of these plots but only to a small degree. Therefore, such edge effects may generally be larger in the treeline ecotone due to the irregularly shaped stems that are seldom vertical. Like the positional errors, edge effects will be relatively larger with decreasing plot size. The effects from trees leaning in and out of the plot, as well as growing inside and outside the plot, will decrease with an increase in the spatial homogeneity of the trees and shrubs.

The graphical display of the AGB se.s (%) at each iteration (Figure 5a–d) indicates that the number of iterations were sufficient for both the TALL and SHORT strata. For all model and stratum combinations, the se.s oscillated around the 0.5% stabilisation criterion after around 10,000 iterations. However, for the SHORT and TALL strata, the se.s values were not stable until after iterations 30,000 and 40,000, respectively. As many as 50,000 bootstrap iterations are quite rare in comparable studies where biomass or other forest attributes are estimated [29,87,88]. The reason is for this is partly because many studies do not use an evaluation criterion for standard error stabilisation. However, the required number of iterations depends on several factors such as the sample size and variability of the original dataset used to fit the model, the complexity of the prediction model used, the desired level of precision of the se estimate, and the uncertainty of the model parameters [89]. In our study, the estimates of the model parameters were characterised by a high level of uncertainty, resulting in a large number of necessary iterations.

Previous studies carried out in productive forests [30,50] have shown that models relying on ALS metrics generally outperform those relying on DAP metrics [52] in terms of prediction accuracy. However, based on the 95% confidence intervals, our study did not find evidence of such a difference. While the R² values indicated weaker relationships between the AGB and DAP metrics compared to the corresponding relationships with the ALS metrics, these weaker relationships were not reflected in the uncertainty estimates. Additionally, the smaller slope parameter values of the DAP models for both strata compared to those of ALS suggest that the DAP point cloud primarily represents the vegetation surface without penetrating the canopies [30,31].

For the modelling, the challenges with the site conditions discussed above meant that we had to choose a robust and simple modelling approach [74]. Even though models for biomass based on point cloud metrics usually include metrics representing both forest height and density, our models only included one statistically significant explanatory variable ($H_{90}$ or $H_{mean})$. Thus, the variation in tree and shrub density, as expressed by the density metrics, was weakly correlated to AGB, and most of the variation could be explained by just the height metrics of our AOI. This might make the models less useful in cases where extrapolation outside the range of variability of the reference data is needed. In our case, however, extrapolation was not extensive as extreme observations were purposely included in our reference data by the way the plots were selected using the ALS data as prior information to guide the sample selection. This was especially important since we relied on linear models [90]. Such a strategy could be adopted in any study area where RS data are available prior to fieldwork.

5. Conclusions

This study indicates that both ALS and DAP data can be used as auxiliary information in area-based AGB estimation with respect to treeline ecotones. We found no significant differences in terms of precision using the two data sources, both at the individual stratum level and for the total area. Although, based on the RMSE and rel.RMSE values, the models for the TALL stratum were more precise than those of the SHORT stratum, the 95% CIs for the models had similar widths in relative terms. This research indicates that area-based biomass estimation can also be carried out operationally in treeline ecotones. However, further studies should be carried out to assess if the proposed method can support repeated surveys for the purpose of AGB change estimation.