Aboveground biomass corresponds strongly with drone-derived canopy height but weakly with greenness (NDVI) in a shrub tundra landscape

We conducted our study on Qikiqtaruk–Herschel Island in the Canadian Arctic. Tundra vegetation communities here range from graminoid- to shrub-dominated and are underlain by organic soils and ice-rich permafrost. This site has undergone marked ecological changes in community composition, increases in canopy height and vegetation abundance, decreases in bare ground, and an advance in leaf emergence and flowering over nearly two decades of ecological monitoring (Myers‐Smith et al 2019).

Our study area encompassed a graminoid-shrub ecotone at the edge of a wet willow shrub-dominated alluvial fan (69.571°N, 138.893°W) (figure 1). We established two adjacent sites, a 0.5 ha area to the north was allocated for monitoring canopy structure and biomass change over time (coordinates in table S4), and an adjacent sampling area to the south contained the harvest plots for calibrating allometric relationships between height or NDVI and biomass. Both the monitoring and harvesting areas were selected to be representative of the full distribution of canopy heights that we observed in the field (see figure S8 for comparison of value distributions) (available online at stacks.iop.org/ERL/15/125004/mmedia).

Zoom In Zoom Out Reset image size

To constrain the photogrammetric modelling and locate the point clouds in a coordinate reference system, 26 ground control markers (265 mm × 265 mm) were deployed across the entire area and geolocated to a relative 3D accuracy of ≤0.015 m with an RTK-GNSS survey instrument (Leica GS10). Coordinates were relative to a local benchmark, geolocated in absolute terms to ±0.003 m in X and Y, and ±0.008 m in Z (95% confidence interval), using the AUSPOS web service. The markers were situated to be visible from the air, and the high density of markers (ca. 26 markers per ha⁻¹) facilitated image alignment for stable reconstruction in the texturally complex scene (Poley and Mcdermid 2020).

In June 2016, we selected 36 plots of 50 cm × 50 cm for detailed observation within our harvesting area (figure 1(d)). The plots were arranged in 12 blocks of three replicates, stratified across the range of canopy heights in order to estimate the allometric models more efficiently as well as to determine the form of the relationship between mean canopy height and biomass (Warton et al 2006). Plots contained no standing water during the period of observations. The dimensions of the harvest plots were selected to be large enough to contain representative samples of plant material and to reduce the possible effects of co-registration errors. The corners of each harvest plot were precisely geolocated using the GNSS. To minimise the GNSS survey staff sinking into the often-soft ground, we used a ca. 25 cm^{2 '}foot' on the bottom of the staff to dissipate pressure.

To enable canopy height to be modelled across the monitoring area, we undertook a walkover survey with the GNSS to measure the ground elevation at intervals along transects, which yielded 911 observations with an average spacing of 1 point per 10 m². We excluded 36 regularly spaced ground observations for validation purposes and interpolated a terrain model from the remaining observations using inverse distance weighting (power = 3, search radius = 7, cell size = 0.1 m).

On the 30th and 31st of July 2016 after the drone surveys were completed, each of the 36, 50 × 50 cm plots were surveyed using point-intercept methods similar to ITEX protocols (Molau and Mølgaard 1996, Myers‐Smith et al 2019). We placed a grid with 36 points at 10 cm intervals over each plot. At each point, we placed a metal pin vertically and recorded the maximum height of the canopy above the moss/litter layer representing the ground surface. Wind speeds were generally light during our ground surveys and our point-intercept observations did not appear to be influenced by the limited movement of the low stature canopies.

Our sampling at the end of July coincided with the period of peak biomass across the growing season at this location. Within each of the 36 sub-plots, all standing vascular plants were harvested down to the top of the moss/litter layer (after Walker et al 2003) on the 31st of July and 1st of August 2016. Harvested biomass was separated into three partitions: (if) woody stems, (ii) shrub leaves (including catkins that accounted for less than 10% of the 'leaf' biomass), and (iii) herbaceous material (consisting of mainly graminoids and equisetum, but also some forbs). Salix richardsonii produces catkins before leaves in June and the seeds are mostly dispersed by mid-July. At the time of drone data collection and biomass harvesting, most catkins had dispersed their seeds and were senesced. Biomass was dried at ca. 35 °C for ≥70 h, until it reached a constant weight (<0.2% change) over a 24 h period.

2.4.1. Aerial survey for canopy height modelling.

2.4.2. Aerial survey for spectral reflectance.

2.5.1. Processing for canopy height models.

2.5.2. Processing for spectral reflectance.

To demonstrate how this approach might support upscaling studies, we estimated aboveground vascular biomass density for the monitoring area based on modelled canopy height and NDVI using the allometric functions calibrated from the adjacent harvest plots (figures 3 and S5, tables S1 and S2). This analysis was undertaken using the 'raster' package v. 3.1–5 (Hijmans et al 2020) and visualised using the package 'rasterVis' v. 0.47 (Lamigueiro and Hijmans 2019). We calculated upper and lower estimates for biomass in the monitoring plot, using the standard error of the allometric equations, in order to account for this source of uncertainty in the landscape estimates.

Statistical analysis was conducted in R (v3.6.1)(R Core Team 2019). To compare agreement between point-intercept and structure-from-motion metrics of canopy height, we calculated concordance correlation coefficients using the 'DescTools' package (after Lin 1989) and we described this relationship with a power function fitted with ordinary least squares regression because using a positive exponent means the model passes through the origin. We used least squares optimisation to fit linear models between canopy height and aboveground biomass, with intercepts constrained through the origin as plants with zero height above ground have no biomass above ground. The intercept terms in the unconstrained point-intercept and photogrammetry canopy height-biomass models were not statistically significantly different from zero (p = 0.78 and p = 0.25, respectively) and constraining model intercepts made only small differences to model slopes (table S1). We reported errors as standard deviations unless otherwise stated.

We used least squares optimisation to fit exponential models between NDVI and three vascular plant biomass pools: (i) total aboveground biomass, (ii) phytomass (calculated as the sum of shrub leaves and herbaceous material), and (iii) the biomass of shrub leaves. Comparisons between remotely-sensed NDVI and biomass usually have a substantial mismatch in observation extents due to the larger grain of satellite products relative to smaller extents of directly measured harvest plots (Berner et al 2018, Bartsch et al 2020). In our study, we undertook spatially explicit drone-based sampling of corresponding areas, so our biomass and NDVI measurements do not have this scale mismatch. Because non-harvested moss can contribute to the differential reflectance of red and NIR energy, we hypothesised that the proportion of moss cover might influence the relationships between NDVI and biomass. We extracted the proportion of moss cover from our point-intercept observations and tested the influence of moss cover on NDVI-biomass relationships by adding an interaction term in our model of the relationship between NDVI and phytomass. The code for statistical analyses and data visualisation is available from https://github.com/AndrewCunliffe/OrcaManuscript.

We found strong agreement between canopy heights as observed with point-intercept method and structure-from-motion photogrammetry (figures 2 and S1). The photogrammetrically derived canopy heights had a consistent positive bias relative to point-intercept heights, with a median difference of 0.14 ± 0.05 m (± SD). Differences in mean canopy height between methods were smaller for the shortest and tallest plots, and greatest for the plots of intermediate heights (figure S1). The concordance correlation coefficient was 0.79 (with 95% confidence intervals of 0.68 to 0.86).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We found canopy height explained most of the variation in the aboveground biomass of vascular plants across the Salix richardsonii-dominated graminoid-shrub ecotone. The biomass models had slopes of 3623 ± 177 g m⁻¹ and 2522 ± 143 g m⁻¹, explaining 0.92 and 0.90 of the variance for point-intercept and SfM-derived canopy heights respectively (figure 3). Total aboveground biomass within the sampled plots ranged from 149 g m² to 2431 g m² with a mean of 1012 ± 699 g m². Shrubs (woody material and leaves) accounted for the majority of biomass in 32 of the 36 plots. The biomass of shrub leaves was positively related to total biomass (slope = 19 g m⁻²) and explained 70% of the variation in total biomass (figure 5(a)). However, phytomass, calculated as the sum of shrub leaves and herbaceous material typically accounted for less than 10% of total biomass (figure S3), and did not correspond with total biomass (figure 5(b)). Herbaceous material (largely equisetum and some forbs) typically accounted for half of the phytomass in each harvest plot, ranging from 3% to 87% of the phytomass. The mass of leaf material was a reasonable predictor of total biomass (y = −63.7 + 19.04x; R² = 0.70; figure 5(a)); however, phytomass was a poor predictor of total biomass (y = −1185−0.8471x; R² = 0.01; figure 5(b)).

We found that NDVI positively corresponded with total aboveground biomass, phytomass and shrub leaf biomass (figures 3(c), 4 and S3, tables 2 and S2). However, NDVI explained less than a quarter of the variance in total aboveground biomass (14% to 23%), phytomass (2% to 7%) and leaf biomass (6% to 21%) across all four spatial grains investigated (figure 4, table 2). The predictive relationships weakened slightly as the spatial grain of the NDVI rasters became finer from 0.121 m to 0.018 m, with larger residual standard errors and smaller coefficients of determination (table 2). With a coarser spatial grain, the overall mean and variability amongst plot NDVI values was lower, although the relationship at the coarsest spatial grain (0.121 m) deviated slightly from this pattern (figure S4(a)). We speculate that this may relate to more pronounced bi-directional reflectance experienced during this particular drone survey that was conducted with a lower sun elevation of just 15.6 degrees (table 1). We tested whether the proportion of moss cover influenced the relationship between NDVI and total biomass, phytomass and the three biomass pools (table S3, figure S6), and though the interaction effects were similar with stronger NDVI-biomass relationships among plots with lower moss cover, the interaction was only significant (p < 0.05) for the phytomass relationship for the 0.121 m raster (figure 5(c)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To understand the performance of NDVI as a predictor of biomass, we examined the pairwise pixel covariance between NDVI and canopy height, as canopy height is a strong predictor of biomass (figure 6). We found a generally positive relationship between NDVI and canopy height, although the relationship was weaker at canopy heights below 0.2 m and once NDVI saturated at ca. 0.75. This relationship was consistent across all four spatial grains of NDVI sampled, and suggests NDVI-biomass transfer functions will be more uncertain where canopies are low (<0.2 m) or NDVI is high (>0.75). We used an NDVI-height transfer function developed by Bartsch et al (2020) to predict canopy height from NDVI and compare these predictions to photogrammetrically derived canopy height (figure S10), and found that while there was a generally positive relationship, the slope differed by a factor of ca. five. Models fitted to predict canopy height as a function of NDVI are reported in table S5.

Zoom In Zoom Out Reset image size

We estimated aboveground vascular biomass density across the graminoid-shrub ecotone for the 0.5 ha⁻¹ monitoring area (figure 1(d)). Estimated biomass at this landscape-level differed substantially between the five predictors, canopy height and NDVI at each of the four spatial grains (figure 6; table 3). We calculated the difference in estimated biomass relative to the CHM-derived biomass map, as we considered this to be the most accurate product (based on the model performance and our knowledge of this site). Although we consider the CHM-derived biomass estimate to be the most accurate, errors in the allometric model and in the canopy heights from the interpolated terrain model and plant surface model will all contribute uncertainty. For example, evaluation of the terrain model accuracy against 36 validation points indicated a mean residual elevation of 0.029 m ± SD 0.035 m. Although small, the systematic bias of 0.029 m in the interpolated terrain model (figure S7) amounts to 10% of the 0.278 m mean canopy height estimated for the monitoring area, and so would equals a 10% underestimate of biomass inferred from the CHM. The allometric models were fitted using similar distributions of canopy height and NDVI values to those found across the adjacent monitoring site. The distribution for the harvested samples of the two coarser NDVI rasters (0.121 m and 0.119 m) had slightly truncated tails, which may have contributed to greater uncertainty in the model fit (figure S8).

We found substantial differences in the estimated average biomass inferred from the five datasets (table 3). NDVI-based biomass estimates are generally greater than the height-based estimates (figures 7(f),(i),(l) and (o); figure S9), apart from areas with taller shrubs near the centre of the monitoring area. Estimated biomass was positively correlated with the spatial grain of the NDVI maps, mainly due to differences in the allometric models between different grain sizes (table S2, figure S5). However, the biomass estimates from NDVI were highly uncertain as demonstrated by the range of possible values (table 3).

Zoom In Zoom Out Reset image size

We found a positive bias in canopy heights measured with point-intercept relative to photogrammetry, which we attribute to differences in the way the two approaches quantify canopy architecture. The photogrammetry-derived heights in our study may have also been slightly exaggerated by slight depression at the plot corners by the survey staff on the ground surface at the top of the moss layer (ca. 2–3 cm) The plot-level terrain models will have smaller vertical errors than the model applied across the monitoring area (figure S7), because elevations were interpolated over smaller horizontal distances (<0.71 m). It is also possible that the fewer point-intercept observations within each harvest plot (n = 36) may under sample canopy structure relative to the photogrammetry (n ≈ 2500). Several studies have now reported strong correspondence between canopy heights reconstructed with photogrammetry and in situ measurements (Alonzo et al 2020, Karl et al 2020, Poley and Mcdermid 2020). For example, Clement and Fraser (2017) reported similarly good correspondence between in situ versus photogrammetrically derived maximum canopy heights for 20 shrubs measured at an Arctic tundra site near Iqaluktuuttiaq (Cambridge Bay). However, such comparisons are hindered by the sensitivity of maximum height measurements to outliers in these often noisy point clouds (Cunliffe et al 2016). The comprehensive review by Poley and Mcdermid (2020) discusses such inter-comparisons of different canopy height measurements in further detail, but they were unable draw general conclusions because all sets of observations are sensitive to the ways in which they are collected, processed and analysed (Cunliffe et al 2016, Wallace et al 2017, Fraser et al 2019). Our findings suggest that, when applied in a consistent manner, drone photogrammetry is an appropriate tool for monitoring shrub canopy heights in such ecosystems.

Canopy height strongly predicted aboveground biomass for the Salix richardsonii-dominated community that we studied, which corroborates similar reports for photogrammetry across a range of biomes and plant communities (Bendig et al 2015, Selsam et al 2017, Grüner et al 2019, Wijesingha et al 2019, Alonzo et al 2020, Howell et al 2020, Karl et al 2020). Estimating aboveground biomass from canopy height models depends on having an underlying terrain model of sufficient quality to describe topographic variability (Cunliffe et al 2016, Fraser et al 2019, Poley and Mcdermid 2020). In this study, we derived our terrain model using RTK-GNSS observations, which can be a viable option for characterising topography over extents of up to a few hectares. In ecosystems where canopies are spatially or temporally discontinuous, terrain models could also be derived directly from photogrammetric point clouds (Cunliffe et al 2016, Fraser et al 2019). Terrain models derived using other survey techniques could also be co-registered in a hybrid approach (Dandois and Ellis 2013). However, propagation of uncertainties including co-registration error is vital for understanding the limits of detection of genuine change in canopy height (James et al 2017).

Relationships between plant dimensions and biomass are sensitive to the ways in which these measurements are obtained (Cunliffe et al 2020). Cross-site data syntheses therefore require the use of standardised protocols for data collection and processing (such as HiLDEN https://arcticdrones.org/, Assmann et al 2018, Cunliffe and Anderson 2019). As noted by Pätzig et al (2020), there is a need for further coordinated work to calibrate the relationship between photogrammetric-inferred canopy height and aboveground biomass for different taxonomic groups. There is also a need to quantify the sensitivity of these relationships to key parameters (e.g. the spatial resolution of the input data, the implementation of multi-view stereopsis and the spatial grain of analysis, sensu Zarco-Tejada et al 2014, Wallace et al 2017), as well as to differences in environmental conditions (e.g. illumination and wind-induced movement of plant canopies, Dandois et al 2015, Frey et al 2018).

We found that NDVI only weakly predicted the aboveground biomass of vascular plants, explaining at most 23% of the variation in total biomass, and even less of the variance in phytomass or leaf biomass (figures 4 and S3, tables 2 and S4). Inferring aboveground biomass from NDVI is predicated on the assumptions that (i) NDVI is a good predictor of phytomass, and (ii) that phytomass is a good predictor of total biomass. We found that while NDVI had some capacity to explain variance in phytomass (figure 4, table 2), phytomass was a very weak predictor of total biomass (figure 5(b)). Across spatial grains, predictive relationships weakened slightly as the spatial grain of the NDVI rasters became finer from 0.121 m to 0.018 m (table 2). We attribute two main causes for the weak correspondence between NDVI and biomass. Firstly, although leaf biomass was a strong predictor of total aboveground biomass, leaf biomass accounted for typically only half of the phytomass in each plot, and phytomass (including herbaceous material and shrub leaves) only weakly corresponded with total biomass (figure 5). Vegetation indices that integrate all photosynthetically active material are often poor predictors of total biomass (Bratsch et al 2017, Räsänen et al 2019). Secondly, we found indications that moss cover influenced the relationship between NDVI and phytomass. The direction of the relationship was consistent; however, the interaction effect was only statistically significant in one of the 12 combinations of NDVI raster and biomass pool tested (table S3). We found that the relationship between NDVI and vascular phytomass was mediated by the amount of moss cover beneath the sampled vegetation and weakened as moss cover increased (figures 5(c) and S6). Vegetation species composition therefore affects the biomass NDVI-relationship.

The relationship between NDVI and biomass in this tundra setting is approximated by the relationship between NDVI and canopy height across the monitoring site. NDVI was generally positivity related to canopy height across all four NDVI rasters (figure 6); however, the relationship was very weak at canopy heights below 0.2 m and NDVI values above ca. 0.75. NDVI-biomass transfer functions will thus be more uncertain where canopies are low (<0.2 m), in some cases due to increased signal from bryophytes, or where NDVI reaches saturation at around 0.75. Canopy heights inferred from NDVI using the NDVI-height transfer function developed by Bartsch et al (2020) were linearly related to canopy heights inferred from photogrammetry, although were positively biased by a factor of ca. five and subject to the same issues around the shorter (<0.2 m) canopies and NDVI saturation. This bias suggests that NDVI-biomass transfer functions will need to be calibrated more specifically for different situations.

The weak correspondence between NDVI and phytomass that we observed contrasts with reports of stronger positive relationships between NDVI and aboveground biomass derived from datasets compiled across different spatial scales (Boelman et al 2003, Walker et al 2003, Goswami et al 2015). NDVI has a saturating relationship with biomass and NDVI-biomass relationships can be confounded by a variety of ecological variables, land-surface properties and view angle effects (Walker et al 2003, Buchhorn et al 2016, Karlsen et al 2018, Myers-Smith et al 2020). Our findings are consistent with the well-known saturation effect (e.g. Berner et al 2018), and we may have found better correspondence between NDVI and biomass if we had sampled over a wider range of NDVI values beyond those that were most prevalent at our study site. The unmeasured biomass associated with moss was likely small compared with the biomass associated with the vascular plants at this site (Reid et al 2012), but if we had included moss in our biomass harvests, this might have modified relationships between NDVI and biomass. Our results highlight the need for caution when deriving total biomass maps from vegetation indices in high latitude ecosystems with variable land cover. The biome-wide tundra greening patterns and trends observed with large-grain satellite datasets are unlikely to directly represent plant functional attributes such as canopy height or biomass in situ (Myers-Smith et al 2020). Thus, to improve our understanding of vegetation greening in tundra ecosystems across vegetation types and geographic gradients, we need data collection across scales from focal sites to the tundra biome (Fisher et al 2018, Miller et al 2019, Myers-Smith et al 2020).

The spatially continuous estimates of canopy height and inferred vascular plant biomass across our monitoring site (figure 6, table 3) are well suited for upscaling studies. These kinds of highly accurate and precise products help to overcome limitations with incomplete characterisations of plant communities due to labour intensive and resource limited in situ monitoring programs (Myers‐Smith et al 2019, Alonzo et al 2020, Bartsch et al 2020). Of particular concern are the difficulties comparing observations from very different spatial extents, between remotely-sensed observations and small, potentially non-representative in situ plots. Biomass estimated from NDVI was highly uncertain (table 3), even when only accounting for uncertainty in the coefficients of our exponential models without accounting for uncertainties in the NDVI values themselves. Drone-derived products are useful for calibrating and validating biomass retrievals of these properties from coarse-scale observations, and for testing key of assumptions underpinning novel retrievals approaches (Bartsch et al 2020). Photogrammetric approaches to monitoring plant canopies can also be deployed over even larger extents using similar data from airborne surveys (Alonzo et al 2020).