Airborne laser scanning reveals increased growth and complexity of boreal forest canopies across a network of ungulate exclosures in Norway

Abstract

Large herbivores are often classed as ecosystem engineers, and when they become scarce or overabundant, this can alter ecosystem states and influence climate forcing potentials. This realization has spurred a call to integrate large herbivores in earth system models. However, we lack a good understanding of their net effects on climate forcing, including carbon and energy exchange. A possible solution to this lies in harmonizing data across the myriad of large herbivore exclosure experiments around the world. This is challenging due to differences in experimental designs and field protocols. We used airborne laser scanning (ALS) to describe the effect of herbivore removal across 43 young boreal forest stands in Norway and found that exclusion caused the canopy height to increase from 1.7 ± 0.2 to 2.5 ± 0.2 m (means ± se), and also causing a marked increase in vertical complexity and above‐ground biomass. We then go on to discuss some of the issues with using ALS; we propose ALS as an approach for studying the effects of multiple large herbivore exclosure experiments simultaneously, and producing area‐based estimates on canopy structure and forest biomass in a cheap, efficient, standardized and reproducible way. We suggest that this is a vital next step towards generating biome‐wide predictions for the effects of large herbivores on forest ecosystem structure which can both inform both local management goals and earth system models.

Introduction

Populations of large herbivores are undergoing drastic changes in many areas, either declining, as is the case for many low latitude areas (Ripple et al., 2015), or increasing, as in many forested ecosystems in the northern hemisphere (Côté et al., 2016; Speed et al., 2019). Consequently, to minimize potential adverse effects of this development, it is paramount to understand the role of large herbivores in their ecosystems. Large herbivore ecologists have long since made this realization (Leopold et al., 1947), and in more recent decades, experimental exclosures have been used extensively for this purpose, as they offer a unique opportunity to study the effects of large herbivore removal on ecosystems (Hester et al., 2000). It is mainly through such manipulative field experiments that we can really get close to and describe cause-effect relationships. By using exclosures, numerous studies have been published on the effects of grazing or browsing animals on diverse aspects of ecology, from biological diversity (Bernes et al., 2018) to ecosystem functioning and climate forcing (Forbes et al., 2019), meaning the potential for herbivores to affect factors that can force climate change, such as carbon dynamics.

Recently, a call has been made to include herbivory in earth system models, as we have realized the potential for large herbivores to influence both carbon storage, sequestration, and albedo (Schmitz et al., 2018). Although work is ongoing, there is still no clear idea about the net effect of large herbivores on climate forcing. For this reason, Leroux et al. (2020) urged the scientific community to come together and create a ‘network of long-term exclosure experiments’, and also point to remote sensing as a tool for gaining a comprehensive understanding of the impact that large herbivores have on carbon source-sink dynamics. We believe these two elements can be united in a joint solution. Analysing data from a diversity of exclosure experiments from around the world is difficult due to varying sampling protocols and experimental designs. This also makes it difficult to find a common denominator between studies which can enable scientists to quantify the climate forcing potential. Remote sensing offers a unique opportunity in this regard, as it can produce area-based statistics in a cheap, efficient, standardized and reproducible way.

Airborne laser scanning (ALS; i.e. airborne light detection and ranging – LiDAR) is an active remote sensing technique which can be used to map both terrain and vegetation in three dimensions. With ALS, laser pulses are emitted from a sensor mounted on an aircraft, and the pulses are distributed perpendicular to the flight line so that a corridor of a certain width is targeted by the emitted pulses. To cover a certain area of interest, multiple parallel flight lines are combined. For each pulse, elapsed time between emission of the laser pulse and return of the backscatter, is used to measure the distance between the sensor and the target (Dubayah & Drake, 2000; Lefsky et al., 2002). In forests, an emitted laser pulse will typically intersect several layers of vegetation, resulting usually in several reflections, or echoes, returning back to the sensor. In multi-layered vegetation, the first echoes represent the upper layer of the canopy. However, since the accumulated backscattered energy needed to trigger an echo has to exceed a certain threshold, the actual point from which the echo is reflected sits below the top of the trees (Thieme et al., 2011). The subsequent echoes from a pulse are reflected from canopy elements closer to the ground, or the ground itself. The point cloud that results from a data acquisition needs to be processed before it can be used efficiently to do analyses on the vegetation structure. First, a digital terrain model (DTM) is constructed. This is enabled by classifying each echo into ‘terrain’ or ‘vegetation’, by means of an algorithm that evaluates the three-dimensional position of the lowest echoes in relation to their neighbouring echoes (Axelsson, 2000). From the echoes classified as ‘terrain’, a triangular irregular network is formed, and the resulting continuous surface serves as DTM, onto which each of the echoes classified as ‘vegetation’ is projected, and their height above the terrain can be derived. The resulting three-dimensional cloud of vegetation echoes provides detailed information on forest structure, such as canopy height, understory foliage density and canopy cover.

Despite its capability to describe vegetation structure in great detail, there are surprisingly few studies that have applied ALS to assess the effects of large herbivores on vegetation (see Davies and Asner (2019) for a notable exception). This could be related to the fact that ALS has been disproportionately used in observational studies (e.g., Melin et al., 2015) over experimental studies, which is how many herbivory studies are conducted. One reason for this must be the mismatch in scale, with ALS data being produced at landscape scales, and experiments being conducted at plot scales (on the scale of meters). However, with ever-increasing resolution (point densities) of ALS data, there is now a possibility of using it at scales down to just a few square meters, making it a viable alternative and cost and time efficient tool compared to traditional field methods (e.g. Thieme et al., 2011).

Norwegian boreal forests are home to some of the world's densest populations of wild cervids, especially the three forest dwelling species moose (Alces alces), red deer (Cervus elaphus) and roe deer (Capreolus capreolus). Although these species are native and with a long history of presence in the region (Grøndahl et al., 2010; Rosvold et al., 2013), this recent development is unprecedented and resulting from a steep increase in numbers post 1970 (Speed et al., 2019). Reasons for this increase in wild cervid numbers include the introduction of directed harvesting, increased use of clear cutting, and reduced competition from livestock (Lavsund et al., 2003). Wild cervids hold a special place in Norwegian culture both as symbols of wild nature and as targets for recreational hunting (Brainerd & Kaltenborn, 2010). Norway is covered by about 37% forest (Statistics Norway, 2021), mainly boreal coniferous forests managed for timber production. Clear-cutting and subsequent planting or natural seeding is the dominant forestry method in this production forest landscape, which appears like a patchwork of different successional stages after logging.

Field experiments using exclosures on recent clear-cuts have revealed a strong effect of moose and other ungulates on the species composition and tree densities in Norwegian boreal forests: across 31 sites, Kolstad et al. (2018a) found that forests protected from ungulates had considerably more large individuals of deciduous species such as birch (Betula pubescens) and rowan (Sorbus aucuparia). However, in this particular experiment, absolute heights were not recorded for trees above 3.5 m and therefore it was not possible to analyse the effect of the fencing treatment on canopy height and other related metrics such as variation in canopy height and vertical canopy diversity. Information on these forest parameters is essential for developing a complete understanding of how herbivores affect climate forcing because of how the canopy height and structure influences energy flows, for example in terms of heat exchange or the absorption of radiation (Thompson et al., 2004). Yet, these parameters are rarely, if ever, reported from ungulate exclosure studies, demonstrating the challenges of developing biome wide experimental networks (as called for by Leroux et al., 2020). Developing protocols for collecting such data is therefore vital.

In this paper we used ALS to describe the effects of ungulate exclusion on canopy height, and variation in canopy height, across 43 boreal production forest sites in Norway, including the 31 sites in Kolstad et al. (2018a). We tested the hypotheses that the cessation of browsing has led to increased forest canopy height and greater canopy structural complexity. We discuss pros and cons of using ALS in experimental ecology as a complement or sometimes as a replacement of more time-consuming field work.

Materials and Methods

This study made use of 45 (later reduced to 43, see below) study sites in three boreal forest regions in Norway (Trøndelag, Telemark, and Akershus/Hedmark; Fig. 1). The dominant tree species in these regions are spruce (Picea abies), pine (Pinus sylvestris), and birch (B. pubescens). Other common trees include rowan (S. aucuparia) and sallow (Salix caprea). Each site consisted of two 20 × 20 m plots which were randomly assigned to one of two treatments: a fenced ungulate exclosure, or an open non-enclosed plot. For sites 17, 18 and 21 (Table 1) the exclosures were 40 × 40 m. The exclosures and the open plots were minimum 20 m apart to avoid edge effects, but still close enough to represent the same climatic and edaphic environments (<100 m). The ungulate fence was 2.5 m tall, denying access to large animals, but allowing small animals like rodents and hares to enter freely. Species of ungulates include wild moose, red deer and roe deer (densities across the study sites are shown in Table 1), as well as free ranging livestock, especially sheep. Livestock, however, only graze during summer; they are less abundant; and they do not usually feed on trees. Therefore, their influence on the tree canopy parameters reported in this study are assumed to be minor. Moose is the most abundant ungulate in terms of metabolic biomass across the field sites, and of the Norwegian cervids it is also the species that depends the most on trees for fodder (Mysterud, 2000). Therefore, it is assumed that differences in tree growth between treatments is primarily, but not solely, a cause of moose exclusion. Fences were erected 0–9 (median 3) years following clear cuttings carried out between the years 2000 and 2009. Early regenerating forests are prime moose habitat, and the reason for selecting these sites was to study the effects of moose browsing on forest ecosystem succession. After initial planting, all silviculture ceased, but forests were still available to the general public for recreational activities.

Figure 1. Map of 45 study sites in Norway. Each point on the map corresponds to the placement of an ungulate exclosure fence and an adjacent open plot. Numbers are site ID’s and refers to Table 1.

Table 1 Site information.

Site ID	Lon.	Lat.	Year when clear cut	Year when exclosures built	ALS from year	Laser point density m−2	Moose density	Red deer density	Roe deer density	Productivity
1	11.3189	64.2637	2006	2008	2010	2	75.89	8.76	7.38	0.22
2	11.8926	64.2134	2005	2008	2011	2	107.57	2.43	17.86	0.12
3	11.8607	64.1997	2004	2008	2011	2	107.57	2.43	17.86	0.09
4	11.3137	64.1813	2004	2008	2010	2	75.89	8.76	7.38	0.08
5	11.5785	63.8439	2006	2008	2015	2	28.27	2.83	6.98	0.11
6	11.6796	63.7363	2003	2008	2015	2	28.27	2.83	6.98	0.24
7	11.7931	63.7219	2005	2008	2015	2	28.27	2.83	6.98	0.01
8	10.7251	63.3583	2002	2008	2016	2	91.60	0.84	20.57	0.16
9	10.4724	63.3522	2002	2008	2017	5	79.89	1.34	31.02	0.20
10	11.0385	63.2549	2004	2008	2015	2	56.45	5.44	4.38	0.27
11	11.0317	63.1931	2003	2008	2015	2	56.45	5.44	4.38	0.12
12	10.9889	63.1360	2002	2008	2015	2	56.45	5.44	4.38	0.12
13	11.4111	63.0981	2005	2008	2015	2	16.48	0.42	0.00	0.21
14	11.4234	63.0470	2005	2008	2015	2	16.48	0.42	0.00	0.13
15	10.9430	62.9673	2005	2008	2015	2	36.26	6.65	2.71	0.04
16	11.9405	60.1946	2008	2011	2019	5	51.74	1.63	10.57	0.51
17	11.5931	60.1692	2004	2007	2019	5	49.09	0.83	16.24	0.50
18	11.7527	60.1245	2004	2007	2019	5	49.09	0.83	16.24	0.27
19	11.9727	60.1229	2008	2011	2019	5	51.74	1.63	10.57	0.10
20	11.8256	60.1093	2004	2007	2019	5	49.09	0.83	16.24	0.16
21	12.1566	60.0956	2007	2010	2016	5	62.65	0.00	13.92	0.34
22	11.6257	60.088	2008	2011	2019	5	52.57	0.36	32.77	0.98
23	11.4894	60.0824	2008	2011	2019	5	52.57	0.36	32.77	0.38
24	11.7035	60.0435	2008	2011	2016	2	52.57	0.36	32.77	0.25
25	11.5109	59.9578	2007	2010	2017	2	56.43	0.00	18.91	1.00
26	11.8852	59.9498	2007	2010	2016	5	62.65	0.00	13.92	0.09
27	11.3105	59.9171	2008	2011	2013	2	59.33	3.44	48.12	0.17
28	11.8016	59.8511	2007	2010	2017	2	56.43	0.00	18.91	0.29
29	11.4273	59.8509	2007	2010	2018	5	56.43	0.00	18.91	0.23
30	11.3548	59.8489	2007	2010	2018	5	56.43	0.00	18.91	0.11
31	11.2610	59.8287	2008	2011	2016	5	71.34	0.00	32.85	0.14
32	8.9705	59.6945	2007	2009	2017	2	45.75	7.55	3.56	0.12
33	9.0807	59.6773	2002	2009	2017	5	45.75	7.55	3.56	0.14
34	9.0119	59.5548	2003	2009	2017	2	45.75	7.55	3.56	0.11
35	9.7385	59.4160	2009	2009	2017	2	45.55	2.23	8.85	0.25
36	8.2668	59.3727	2000	2009	2017	5	34.25	13.40	2.94	0.09
37	9.8235	59.3672	2005	2009	2017	2	45.49	33.90	19.92	0.15
38	8.1985	59.3433	2004	2009	2017	5	34.25	13.40	2.94	0.11
39	8.0925	59.3370	2007	2009	2017	2	23.22	02.08	1.53	0.06
40	8.5493	59.3092	2005	2009	2017	5	34.25	13.40	2.94	0.11
41	9.1249	59.2619	2006	2009	2017	2	30.19	28.26	15.05	0.17
42	9.2637	59.1417	2006	2009	2017	2	36.73	7.24	26.88	0.08
43	9.1470	59.1279	2005	2009	2017	2	27.01	9.84	05.03	0.21
44	8.8338	59.0570	2007	2009	2016	2	27.01	9.84	05.03	0.02
45	8.8225	59.0217	2007	2009	2016	2	27.01	9.84	05.03	0.21

Site ID refers to numbers in Figure 1. Cervid densities are metabolic biomass (kg km⁻²) estimated at the municipality level for year 2015, reported in Speed et al. (2019). Productivity is a scaled variable indicating tree production potential, and is calculated as in Kolstad et al. (2018a). Shaded rows represent sites that were removed from analyses as they were seen as outliers along the productivity gradient.

The study sites were originally selected with the intent of covering large variation in forest productivity, and in (Kolstad et al., 2018a) we developed an index of productivity which we here use as a covariate in our statistical models. This index is a scaled (between 0 and 1) and continuous variable reflecting the mean annual tree biomass increment, which we obtained by following individual trees and measuring their height every year, and converting this to estimates of biomass using allometric equations. The productivity index was developed initially for the Trøndelag region (Kolstad et al., 2018a), but was updated to include the sites in the Akershus and Hedmark region in the current study (Fig. 1; Table 1). Two of the sites (number 22 and 25 in Table 1) were removed from our subsequent analyses as they were considered to be extreme with regard to the productivity variable (2.45 and 2.37 sd-units from the third highest value of productivity).

To establish a new, detailed terrain model for Norway for topographic mapping and management that requires detailed terrain information, ALS data for the entire land area below the elevational treeline will be freely available by the end of 2022. The program is called ‘Nasjonal detaljert høydemodell’, NDH). For our study sites, data were already available and downloaded from hoydedata.no (Kartverket, 2019) as a 1 × 1 km square around each study site. We selected the most recent data if there was repeated sampling, which gave an experimental duration (years between when fences were erected and when ALS data were collected) of between 2 and 12 years (median 8). Point densities were either 2 or 5 m^-2, which have proved to be sufficient for forest inventories of mature stands (Gobakken & Næsset, 2008), and for the prediction of height and stem densities in regeneration forest (Ørka et al., 2016b). The ALS data were recorded using pulse sensors with full waveform recording; the most used sensor in NDH is Riegl Q1560i. The laser pulse footprint is generally smaller than 0.61 m², depending on flight height. The data collection platform was either a fixed wing aircraft or a helicopter.

The data was decompressed from LAZ to LAS format and imported into R (R Core Team, 2019) [version 3.6.1.] where all subsequent analyses were carried out, using many of the functions from the lidR package (Roussel & Auty, 2020). The LAS-files (ALS point clouds) were initially clipped to include a 6 m outer buffer around each plot to reduce edge effects of the subsequent individual tree detection (see below). When clipping the data extent we made use of plot coordinates recorded by standard handheld GPS with an accuracy of c. 6 m. To correct initial inaccuracies in plot coordinates, the coordinates for the exclosures were moved so that they overlaid the corners of the fence observed in aerial orthophotos. Since there was no marking of the open plots visible from the aerial photos, the coordinates of the open plots were moved an equivalent distance and direction as the coordinates of the exclosures. We believe this to be justified, as the errors for the coordinates in the exclosures appeared on visual inspection to have some systematic component, but we acknowledge that this is a source of uncertainty. Note also that the open plots were generally situated in homogeneous clear-cut stands, reducing the potential consequences of georeferencing errors.

Some larger trees were left standing after the clear cutting and these had to be removed from the point cloud before any further analyses because they were not accessible for cervids to browse on and can therefore be seen as not part of the experiment. To do this, first, the height of the ALS points were normalized to represent height above the ground (rather than above sea level) by creating a DTM and subtracting the height values from that, using the normalize_height function in the lidR package. Second, the find_trees function was used, with a moving square window of 8 m, to identify individual trees. The points associated with each tree was classified using the segment_trees function [dalponte2016 algorithm (Dalponte & Coomes, 2016): growing threshold 1 = 0.35; growing threshold 2 = 0.55; maximum crown diameter = 6 m]. This function also makes use of a canopy model (raster file), which was produced using the grid_canopy function and a point-to-raster method.

The LAS files were clipped, using the clip_roi function, to remove trees taller than 10 m (the remnant trees left standing after clear cutting), and subsequently to remove the outer buffer around the experimental plots, and another 2 m inner buffer, resulting in a plot size of 18 × 18 m. Removing data 2 m in from the plot edges was done to minimize the chances of ALS data not overlapping with the experimental plots, and also to minimize ecological edge effects, such as altered light and wind penetration (Fig. 2). The final point clouds after filtering contained between 405 and 9486 vegetation echoes (mean=1883).

Figure 2. A canopy height model for the open plot at site 9 showing the extent of the spatial data which is 32 × 32 m, the exclusion fence (20 × 20 m, corners marked with circles) and the final extent after removing a 2 m inner buffer (18 × 18 m). The green areas are remnant tree left after clear cutting, and these were removed from the ALS point cloud by cutting around the outlining polygons.

Multiple metrics representing both the vertical and horizontal distribution of the point clouds were computed using the first echoes only (Ørka, 2016a). To represent canopy height, the 95th percentile was used (Fig. 3). This was instead of the maximum value as that is more sensitive to variation in laser point densities (see Gobakken & Næsset, 2008). As a measure of vegetation heterogeneity we calculated the vertical complexity index (VCI). The VCI is a scale-independent measure of evenness and is 1 for a complete uniform distribution, and decreases as the point cloud becomes more skewed (van Ewijk et al., 2011). Above ground biomass (AGB; Mg ha⁻¹) values for our 18 × 18 m plots (324 m²) were predicted using the equation for AGB1999 in Økseter et al. (2015), dependent on ALS metrics. The equation as it is written in Økseter et al. (2015) contained a typographical error which we accounted for (intercept should read ln(35.55); Økseter 2021, pers. comm.). This model was constructed using field plot (200 m²) observations from young coniferous stands aged between 8 and 23 years. There are several concerns related to the application of external models dependent on ALS metrics. Studies have shown that both variable selection and the values of model parameter estimates are dependent on both the properties of the field data and technical issues such as the laser instrument used, and the setup used for the data acquisition (de Lera Garrido et al., 2020; Næsset, 2005). Furthermore, ALS metrics are scale dependent (Packalen et al., 2019) which means that it is not ideal to apply models constructed for basic units of 200 m² to units of 324 m². However, for the purpose of this study, the possible shortcomings will most likely not affect the conclusions. We acknowledge that there is a risk that our predictions can have systematic errors, but we expect the level to be similar for each of the sites and consistent between sites.

Figure 3. ALS canopy height distribution for an 18 × 18 m ungulate exclosure and an adjacent open area. The data comes from one of the most productive sites (site 17), where the experimental duration has been 12 years. Arrows indicate the position of the 95th percentile which is used as a measure for canopy height in this study. The figure illustrates how the foliar density (or density of points) inside large herbivore exclosures has increased in height above the ground.

Linear mixed effects models (Brooks et al., 2017) were used to investigate the effect of ungulate exclusion and its interactive effects with site productivity, and experimental duration, on three response variables: canopy height, VCI, and AGB. A random intercept was added for each site. Random intercepts for the three regions were found to not improve the model fit and were therefore not included. The response variables and the continuous predictors were divided by 1 standard deviation (sd) and centred on the mean value, whereas the binary predictor ‘Treatment’ was coded as either −0.5 or 0.5, and thus with a mean of zero. Global models also included a quadratic term for site productivity, as well as the interaction terms treatment × experimental duration, and treatment × productivity. The global models were run through the dredge function in the MuMIn package (Bartoń, 2020) to find all possible model configurations and rank them according to a criterion of lowest AIC_c score. Models that scored less than 2 AICc units from the best model were included in a confidence set of models, except for models that included the quadratic term for productivity but not the linear term. Model averaging (Bartoń, 2020) was done for each confidence set based on the proportional variance explained by each model, using the subset, or conditional, method (Grueber et al., 2011). Models were visually validated for all model assumptions.

Results

The canopy height (the 95th percentile) varied between almost zero to around 6 m; it was higher inside exclosures and this difference increased with time since exclusion (Figs. 4A, B and 5). The mean (±se) canopy height (m) was 1.73 ± 0.2 and 2.46 ± 0.2 for the open plots and the exclosures, respectively, corresponding to a 42.7% increase following ungulate exclusion. The relative variable importance for the main effect of herbivore exclusion, and its interaction with experimental duration, were both 1 (terms appearing in all models of the candidate set with deltaAICc <2.0), and the effect sizes were 0.51 ± 0.12 and 0.37 ± 0.12 sd-units, respectively (Fig. 5; Table S1). Productivity was not an important predictor for canopy height, neither by itself, or through an interaction with herbivore exclusion (Fig. 5).

Figure 4. Variables of interest derived from airborne laser scanning from 43 boreal forest sites in Norway which each had one 20 × 20 m ungulate exclosure and a corresponding and adjacent open plot. (A) Canopy height (defined as the 95th percentile) increased with experimental duration, and this increase was greater inside exclosures. Dotplots showing the distribution of canopy height (B) and above ground biomass (C) across the 43 study sites and the two treatments. The diamond symbol represents the mean. (D) Vertical complexity index (VCI) is higher inside exclosures and shows a diminishing positive relationship with site productivity. Regression lines are predicted values (±se) derived from averaged linear mixed effects models, with all other explanatory variables set to their mean value. Note: AGB is predicted from an external model and hence the units on the y-axis are therefore approximate.

Figure 5. Results from multi-model inference on three response variables: above ground biomass (AGB), Vertical Complexity Index (VCI) and Canopy Height. Relative variable importance (horizontal bars; left) and standardized, model averaged parameter estimates with 95% confidence intervals (right) are shown. Predictors that include ungulate exclusion are of special interest and are coloured red. The response variables are also scaled, allowing a direct comparison of effect sizes between the two models. The predictors are centered, including the categorical predictor Herbivore Exclusion, so that main effects are interpretable as their effect when all other variables are set to their mean value, and irrespective of herbivore exclusion. The numbers behind the figure are also presented in Table S1.

AGB ranged from almost zero to 39.4 Mg ha⁻¹, with a median of 2.4 and mean 4.9; it was higher inside exclosures and this difference increased with time since exclusion (importances = 1; Figs. 4C and 5). The mean AGB (Mg ha⁻¹) was 2.91 ± 0.7 and 6.83 ± 1.2 for the open plots and the exclosures, respectively, corresponding to a 3.9 Mg ha⁻¹, or 235%, increase following ungulate exclusion.

There was a positive effect of ungulate exclusion [relative variable importance = 1; effect size = 0.6 (± 0.1) sd-units], productivity, and experimental duration on the VCI, but no strong interaction effects between ungulate exclusion and productivity or experimental duration (Figs. 4D and 5).

Discussion

In this study we used ALS to study the effects of 43 ungulate exclosures on boreal forest canopies across Norway. Our results show that wild cervids have a considerable effect in reducing forest canopy growth following an initial disturbance. Our study also shows that excluding ungulates increased both the vertical canopy diversity and the AGB (hence also aboveground carbon stocks) of the early regrowth forests. This study provides support for using ALS to analyse diverse ungulate exclosure experiments across biomes (see Vuorinen et al., 2020) in a standardized and cost-effective way, and puts us closer of being able to quantify the role of large herbivores in the global climate system (Leroux et al., 2020; Schmitz et al., 2018).

The canopy was 43% taller inside exclosures compared to the open plots. This is in support of our hypothesis, and is in line with previous analyses of the same study sites (Kolstad et al., 2018a), and with other ungulate exclosure studies in boreal forests (Conway & Johnstone, 2017; Hidding et al., 2013; McInnes et al., 1992). Similarly, the AGB was 235% greater inside the exclosures, which corroborates the field based results from a subset of the same ungulate exclosures where the above ground carbons stock was found to be greater inside exclosures (Gebreyohanis, 2017). Note, however, that because these stands are very young, the absolute difference between treatments is still moderate (3.9 Mg ha⁻¹), meaning that the significance of this difference in terms of its global climate forcing potential depends on whether the relative difference in AGB between treatments is maintained as the stands develop, for example due to long-lasting changes in species composition (Hidding et al., 2013). It is therefore of great interest to be able to extrapolate our findings into the future. However, these findings are somewhat context dependent and are not easily transferable to older forest stands, nor to forests subjected to commercial thinning. This is primarily because of two factors:

Firstly, moose, the dominant herbivore in these regions, feed primarily in early successional forests (Bjørneraas et al., 2011; Wam et al., 2016), and less so in mature forests. Therefore, the exclusion effect can be expected to decrease over time, possibly becoming reversed, if the open plots close up with unpalatable spruce (Månsson et al., 2007), which will grow taller than most other species in these forests. Rather than impeding forest growth as such, it is perhaps more correct to interpret this effect as the potential for wild cervids, especially moose, to suppress the initial deciduous stage following clear cutting (Kolstad et al., 2018a; Vuorinen et al., 2020). This is in accordance with existing theory that browsing in unproductive ecosystems (such as boreal forests) will cause an acceleration of succession towards late successional species (Davidson, 1993), in this case spruce.

Secondly, as our experimental plots were not subject to pre-commercial thinning (also called cleaning), fertilization, or pesticide treatments, generalizations to commercially driven production forests can be difficult. Therefore, in order to fully explore the effects of large herbivores on forest development, exclusion studies should preferably continue for the entire life cycle of a forest stand and include interactions with the forestry practiced in the area. Subsequently, or with these issues otherwise addressed, for example by creating mechanistic models linking above ground characteristics to soil carbon storage (Kolstad et al., 2018b), one can attempt to scale up the effects of browsing by wild cervids to much larger extents, such as countries or biomes (e.g. Schmitz et al., 2014).

We predicted AGB using an external model constructed with data from an ALS acquisition outside our study area. This introduces the possibility of systematic errors on the predicted values (de Lera Garrido et al 2020; Packalen et al., 2019). We considered the relative treatment difference to be little affected by these potential systematic errors caused by the model, although an aggregation of predictions for the different sites would be biased and the errors are unknown. The absolute values of AGB (reported as Mg ha⁻¹) should therefore be treated as approximate. To overcome this, if the goal is to estimate AGB for an area, we recommend using field data to calibrate the model within each ALS acquisition (de Lera Garrido et al., 2020).

Excluding ungulates increased the vertical canopy diversity across our study sites, indicating that the canopies inside the exclosures varied more in height compared to those in the open plots. This matched our hypothesis, and corresponds well with the observation that the palatable deciduous trees in the open plots can be so heavily browsed that their heights converge around the approximate average snow depth, forming a type of browsing lawn (Fornara & Du Toit, 2007). Our metric for analysing this type of structural canopy complexity was the VCI, but there are many more to choose from (e.g. van Ewijk et al., 2011), each with its own merits. Increased vertical canopy diversity will have implications for heat and gas exchange, leading possibly to localized warming due to increased sensible heating (Thompson et al., 2004). However, changes in forest structural complexity is a rarely studied consequence of ungulate herbivory (but see Eichhorn et al., 2017). Increased structural diversity may also influence biodiversity and habitat suitability for forest dwelling species, such as birds (Cardinal et al., 2012; MacArthur & MacArthur, 1961), but these relationships are predominantly studied in mature forests, and more studies from early successional forests are needed.

Although this study focuses on canopy aspects relevant for understanding the climate forcing effect of browsing ungulates, the by far largest carbon pool in this ecosystem is in the soil, and that is not directly detectable by remote sensing. In a study based on the 15 sites in the Trøndelag region 8 years after exclusion of ungulates, Kolstad et al. (2018b) did not find a treatment effect on soil carbon stocks, which were on average 87.17 Mg ha⁻¹ (compared to 2.4 Mg ha⁻¹ above ground carbon estimated in this study, assuming a 50% carbon content of wood). Remote sensing is therefore not a singular solution for evaluating the effect of large herbivores on ecosystem carbon stocks.

Large herbivore exclosure experiments are generally not designed with the intention of using remotely sensed data as auxiliary information. For example, precise GPS coordinates are perhaps not considered essential, as the exclosures are easily seen in the field and open plots are usually marked with stakes. This will, however, introduce georeferencing errors. Also, large herbivore exclosures are usually quadratic, as in our case, or sometimes rectangular. If relying on coupling field data with remote sensing data, any sampling shape differing from a perfect circle will increase the perimeter-to-area ratio and therefore also increase the possibility for co-measurement errors. Adding an internal buffer of a few meters can alleviate this to some extent, as we did in this study.

The median experimental duration of this study was 8 years; although the exclosures are considerably older in some cases, the experimental duration is contingent on the year of the ALS acquisition. Because coordinated ALS campaigns are not conducted frequently, very recent data cannot be expected, and rarely can baseline data, that is, data from before the experiment started, be obtained. This is a major constraint on the use of ALS for experimental field plots. Using terrestrial laser scanning (e.g. Eichhorn et al., 2017), or an ALS system attached to a drone, can overcome this problem, but both these platforms have limitations with regard to spatial extent. However, with increasing length of the field experiment, more and more ALS data will become overlapping with the field study, allowing the analyses of time series. As an alternative to ALS, three-dimensional point clouds can also be obtained from aerial images using digital photogrammetry (White et al., 2013). The collection of aerial images are more frequent than that of ALS, and since such point clouds have shown to accurately depict the canopy surface when a DTM is already available (Gobakken et al., 2015), aerial images could possibly be a great source of large-scale and temporally resolute remotely sensed data.

In conclusion, excluding ungulates from early successional boreal forests in Norway led to a marked increase in canopy height, vertical canopy complexity and AGB. This indicates some of the potential for wild ungulates to significantly influence the local successional pathways and simultaneously the global energy balance. Furthermore, this work demonstrates how ALS could be used for understanding the ecological effects of large herbivores on forested ecosystems, with implications for climate forcing and components of biodiversity, and how it could potentially be used to analyse treatment effects from diverse large herbivore exclosure experiments, regardless of variations in field sampling protocols. We see this as a vital next step towards generating biome-wide predictions for the effects of large herbivores (as called for by Leroux et al 2020) which can inform both local management goals and earth system models.

Acknowledgements

We thank Marc Daverdin at the NTNU University museum for preparing the raw LAZ files for us. We are also grateful to the landowners who have allowed us to establish this long-term experiment in their forests. The Sustherb and Disentangle projects received funding from the Research Council of Norway (Project 184036 & 262064), the Norwegian Environment Agency, and the County Administrations and County Governors in Trøndelag, Møre og Romsdal, Vestland, Innlandet, Viken and Vestfold og Telemark counties. Finally, we wish to thank Shawn J. Leroux and one other anonymous reviewer who critically assessed the manuscripts and gave constructive feedback.

Funding Information

The Sustherb and Disentangle projects received funding from the Research Council of Norway (Project 184036 & 262064), the Norwegian Environment Agency, and the County Administrations and County Governors in Trøndelag, Møre og Romsdal, Vestland, Innlandet, Viken and Vestfold og Telemark counties.

Data Availability Statement

The dataset for the analyses presented here are available online (Kolstad et al., 2021).