Insights into Boreal Forest Disturbance from Canopy Stability Index

Abstract

The world’s forests are being increasingly disturbed from exposure to the compounding impacts of land use and climate change, in addition to natural disturbance regimes. Boreal forests have a lower level of deforestation compared to tropical forests, and while they have higher levels of natural disturbances, the accumulated impact of forest management for commodity production coupled with worsening fire weather conditions and other climate-related stressors is resulting in ecosystem degradation and loss of biodiversity. We used satellite-based time-series analysis of two canopy indices—canopy photosynthesis and canopy water stress—to calculate an index that maps the relative stability of forest canopies in the Canadian provinces of Ontario and Quebec. By drawing upon available spatial time-series data on logging, wildfire, and insect infestation impacts, we were able to attribute the causal determinants of areas identified as having unstable forest canopy. The slope of the two indices that comprise the stability index also provided information as to where the forest is recovering from prior disturbances. The stability analyses and associated spatial datasets are available in an interactive web-based mapping app. that can be used to map disturbed forest canopies and the attribution of disturbances to human or natural causes. This information can assist decision-makers in identifying areas that are potentially ecologically degraded and in need of restoration and those stable areas that are a priority for protection.

1. Introduction

Forests yield a range of ecosystem services that provide significant social, environmental, and economic benefits for people, including climate mitigation, clean water, and wildlife habitat [1]. However, the world’s forests are being increasingly disturbed from exposure to the compounding impacts of land use and climate change, in addition to natural disturbance regimes [2,3,4,5,6]. Boreal forests have a lower level of deforestation compared to tropical forests. While they experience higher levels of natural disturbances, the accumulated impact of forest management for commodity production—coupled with worsening fire weather conditions and other climate-related stressors—has been shown to result in various forms of ecosystem degradation. This includes failed regeneration [7], loss of biodiversity [8], loss of stand age diversity, particularly older forest, and loss of critical caribou habitat [3]. However, a disturbed canopy does not necessarily equate with a particular kind of degradation or a degraded forest per se, as the ecological consequences depend on the type, intensity, frequency, duration, and geographic extent of the disturbances, i.e., the disturbance regime [9]. To understand if and how forests are being degraded, information is first needed on how boreal forests are being disturbed, whether the causes are natural, anthropogenic, or some combination of factors, as well as their accumulated impacts geographically and over time. This information provides important context for further monitoring and evaluation of the ecological condition of forests that can be used to guide forest planning and management responses [10].

It is also important that information on forest stability and the attribution of unstable sites to human or natural causes is available to conservation practitioners, including non-government and government organizations. A recent survey and workshop identified information needs to move from knowledge to action for biodiversity conservation in Canada [11]. The findings focused on understanding major threats and improving mobilization accessibility of information through: an improved understanding of approaches to monitor and mitigate the cumulative and interactive effects of multiple stressors that alter the ability of habitats to support biodiversity; development of methods to integrate data across regions, taxa, and spatial-temporal scales; and development of approaches to improve access for practitioners, managers, and the public to remote sensing and GIS biodiversity data of consistent quality across sources and disseminate resulting evidence. Some of these information needs can be met through the development of an interactive web-based mapping tool that provides users with a low level of spatial analytical technical skills and the ability to undertake exploratory space-time analysis of forest canopy stability and where unstable canopies are associated with natural and/or anthropogenic disturbances. Determining the underlying factors that drive forest stability is important for science and conservation [12].

Advances in satellite, airborne, and drone remote sensing, coupled with machine learning trained on networks of field observations, are being exploited to generate a growing array of Earth system data. These data products are used to monitor changes in forest cover and structure at increasingly fine spatial and temporal resolutions [13,14,15,16,17]. Furthermore, they are becoming more accessible to a wider range of users through web-based interactive mapping tools [WRI-GFW reference]. While valuable for many real-world applications, including mapping forest loss, these global forest-related modeled data have their limitations, including representing only a “snapshot” in time, having spatially varying accuracy due to limited field observation data for model calibration, and being unable to attribute a disturbance to anthropic or natural causes. Some of the limitations can be addressed by using indices based on time-series analysis of satellite data which can capture temporal gradients in forest conditions such as forest canopy stability indices [18].

Here we explore the question of whether the limitations of Earth systems data for assessing forest disturbance regimes, attributing causal factors, and providing context for further studies on forest degradation, can be addressed with complementary data sources and specifically (1) time-series analyses that can identify past disturbances and whether a forest area is in recovery and (2) regionally sourced land use, wildfire, and insect infestation records. We proposed that these complementary data sources provide insights into the current ecological condition of a forest and the attribution of disturbances to anthropogenic or natural causes.

Our investigations were focused on the boreal forests of two Canadian provinces (Ontario and Quebec) as a case study. We first derived an index of forest canopy stability from a time-series analysis of Earth system satellite-based data, which we used to map disturbed areas. We then accessed and compiled spatial databases of long-term (~50 years) forest management, wildfires, and insect infestations. We then undertook an exploratory analysis using machine learning to identify the causal factors that are correlated with areas identified as having experienced disturbances.

2. Materials and Methods

We followed a structured workflow to derive an index of forest canopy stability and explore the disturbance factors, as per the overview provided in Figure 1. This involved two major steps: (i) calculating the forest canopy stability index based on time series of remotely sensed canopy indices, and (ii) identifying the disturbance factors causing unstable forest canopy (i.e., low stability index values). The workflow was applied to all forests and just to boreal forests [19] within the Canadian provinces of Ontario and Quebec, giving a total study region of around 50 M ha.

2.1. Step 1—Forest Canopy Stability Index

Further details on the analytical steps involved in the first step—calculating the stability index—are shown in Figure 2. First, two canopy-relevant indices are calculated at a 500 m pixel resolution from the MODIS satellite-based sensor for growing seasons months for the years 2003–2022. From the monthly values, the Coefficient of Variation and the slope of change for each of the two indices give four index values for each pixel, which are ranked into deciles, summed, and normalized to produce stability index values of 1–10, with 10 being the most stable forest canopy. Further details on each step are provided below.

The forest canopy stability index was calculated using a method developed by Shestakova et al. [18], building on the approach of Mackey et al. [20]. We created a time series from 2003 to 2022 for two metrics correlated with canopy photosynthesis rates and desiccation: (1) the fraction of photosynthetically active radiation intercepted by the canopy (fPAR), and (2) the shortwave infrared water stress index (SIWSI). These metrics were derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) [21,22] on a mean-monthly basis with a 500 m pixel resolution, analyzing only good quality pixels (quality flag 0) with clear cloud states and excluding pixels with snow/ice, cirrus, internal cloud masks, and cloud shadows.

To account for different forest types and phenology, we spatially stratified this study region, focusing solely on forested areas. Non-forest areas were excluded using a 30 m pixel resolution land cover dataset from 2019 [23]. We further stratified forest cover into coniferous, broadleaf, and mixed-wood areas, downscaling and reprojecting the data to match the MODIS projection in Google Earth Engine [24] (Figure 3a). The command ‘reduceResolution’ was used to downscale the data with the mode reducer and reproject using the ‘reproject’ command. The stability index calculations were confined to analysis of fPAR and SWISI during the months of the year that correspond to the growing season. However, the start, end, and length of the growing season vary across this study region with latitude. Thus, we classified pixels by their growing season start and end weeks based on 2020 gridded data [25] (Figure 3b).

Modifying the index of Shestakova et al. [18], we calculated the stability index using four canopy metrics based on monthly mean values from 2003 to 2022: (1) Coefficient of Variation (CoV) of annual mean fPAR; (2) absolute slope of annual mean fPAR; (3) CoV of annual mean SIWSI; and (4) absolute slope of annual mean SIWSI. Each metric was binned into percentiles, assigned class values from 1 to 10, and summed to produce class values ranging from 4 (most stable) to 40 (least stable). These class values were then divided by 4 and subtracted from 11 to give index values ranging from 1 (least stable) to 10 (most stable). The index was calculated for each of the 138 spatial units (derived from the unique combinations of three forest types and 117 growing season values) with the same forest cover type and growing season region, then mosaicked to create a stability index map for forests in Ontario and Quebec (Figure 4). All calculations were performed in Google Earth Engine.

2.2. Step 2—Identifying Disturbance Factors That Cause Forest Canopy Instability

The second step in the analytical workflow involved the following: (i) assembling datasets on disturbance factors; and (ii) conducting Gradient Boosted Machine (GBM) modeling to identify the disturbance factors associated with unstable forest canopies (i.e., pixels with low stability index values); and examining two diagnostic statistics (practical dependence plots, PDP; individual conditional expectation plots, ICE) to identify the optimum model for current purposes.

2.2.1. Disturbance Factors

We identified a set of variables that could serve as potential explanatory drivers of forest instability, for which there were available spatial data, and which were either direct forest disturbances, associated with disturbances, or known to influence the intensity of disturbances:

(a)

Direct drives of forest disturbance (harvest/logging, wildfire, insect infestation)—Data on the presence/absence of wildfire from years 1973 to 2020 [26,27,28,29,30] and forest harvesting from years approximately 1976 to 2020 [3,26,30,31,32] were compiled from national and provincial sources. Insect disturbance severity was ranked and combined from provincial datasets (Figure S1) [33,34]. Insect disturbance layers were combined and ranked for Ontario and Quebec as shown in Tables S2 and S3.

(b)

Indirect correlates of forest disturbance (land use and infrastructure; roads)—Land use and infrastructure disturbances included powerlines, railways, seismic lines, pipelines, dams, airstrips, mines, reservoirs, settlements, well sites, agriculture, and oil and gas, derived from provincial layers [3,35,36]. The road network [37,38,39] was buffered to 500 m and included as a separate input. Numerous studies have shown that proximity to roads and urban settlements, along with other built infrastructure, impacts forest cover, structure, and composition and elevates the probability of fragmentation [40,41,42,43].

(c)

Environmental factors associated with different disturbance regimes (landform classes, geological substrate, Topographic Wetness Index (TWI)). Previous studies have shown that these environmental factors influence the likelihood or intensity of wildfires, harvesting, or insect infestations [44,45,46,47]. The Global SRTM Landforms [48] were simplified into categories as shown in Table S1. Surficial geology was categorized based on Canadian Geoscience Maps [49], as seen in Table S4. TWI describes the topographic controls on the tendency for an area to accumulate water and was calculated [50] using the SRTM digital elevation model [51] and Whitebox tools [52].

2.2.2. Building the GBM Models

To investigate which disturbance factors are most strongly associated with unstable forest areas (i.e., pixels with a low stability index), we used the GBM algorithm, employing the gbm function from the gbm package [53] in R [54]. The GBM algorithm, which is suitable for modeling complex non-linear relationships, built an ensemble of decision trees to predict forest instability. The dependent variable was the forest stability index, with the potential explanatory variables being fire; harvest; insect infestation; landform; (non-forestry) land use; roads; geology; and TWI. Over one million random points (20% of the pixels in our study area) were generated using ArcGIS Pro [55] to sample these datasets for input into the GBM algorithm in GRASS GIS [56]. All layers were downscaled, reprojected, and clipped to match the stability index layer using ArcGIS Pro. Using the “Resample” tool in the Data Management Toolbox, the ‘majority’ algorithm was used to downscale to the stability index resolution and was reprojected in the environmental variables. We split the data into training (70%) and testing (30%) datasets, tuning the GBM model to find optimal hyperparameters: interaction.depth (1, 3, 5), shrinkage (0.01, 0.05, 0.1), bag.fraction (0.5, 0.75, 1), and n.minobsinnode (5, 10, 15). Using 5-fold cross-validation, we evaluated 81 hyperparameter combinations, selecting the model with the lowest root mean squared error (RMSE). There were four model runs: (i) all forests and all explanatory variables, (ii) all forests and just the environmental variables, (iii) boreal forest and all explanatory variables, and (iv) boreal forest and just environmental variables. The final model parameters for each of the four runs can be seen in Table S5. The relative influence values of the independent variables are outputs of the chosen model calculation.

2.2.3. Partial Dependence Plots (PDP) and Individual Conditional Expectation (ICE) Plots

PDP and ICE plots were used to visualize the relationship between each variable and forest stability. The PDPs show how the stability of the forest changes as the variable of interest changes while keeping all other variables constant, indicating the importance of variable of interest [57]. PDP shows the average of the effect of the variable on stability, while ICE shows the relationship between stability and variable of interest for an individual sample. ICE plots highlight the variation in the distribution of the variable effect on stability [58]. For the continuous variable, TWI, the ICE plot was centered to a first PDP value and displays the differences in prediction for easier interpretation. The PDP and ICE values were calculated using the ‘partial’ function from the pdp package [59] in R.

2.2.4. Ablation Study

An ablation study was conducted to assess the impact of removing each variable from the model. We reran the model, excluding one variable at a time, and recorded the RMSE, optimal number of trees, mean absolute error (MAE), and RMSE of the predictions. This analysis, along with the GBM modeling, helped determine the importance of each variable in predicting forest instability. The GBM modeling and ablation study were performed on the following: (a) boreal forests as a function of (i) all potential explanatory variables and (ii) a subset of environmental explanatory variables; and (b) for all forests in this study area (i.e., the boreal forest plus the non-boreal forest to the south) as a function of (i) the full set of independent variables and (ii) a subset of environmental variables.

2.3. Time Series Pixel Drill

We built an interactive, web-based mapping interrogation tool that enables pixel drills at any user-specified location within this study area. For each selected pixel, a time-series graph is plotted for the fPAR and SIWSI, the two-component indices of the stability index, along with data on known logging, wildfire, and insect infestation events. The tool can be accessed at the website https://www.globalforestobservatory.com/ (accessed on 29 September 2024).

3. Results

The full range in the stability index is evident across this study region (Figure 4). Extensive areas of highly unstable forest are evident in the far west and Southeast, while a scattering of smaller areas of instability is evident. Areas of high forest stability are found in the north and around the center (Figure 4).

In terms of causation, we focus here on the results of the GBM model run that was focused on boreal forests using the full set of potential explanatory variables. The results for the “all forests” model (which included boreal plus the non-boreal forest to the south) and for the two model runs based on just the environmental variables can be found in Supplementary Materials. The relative influence of each potential explanatory variable on the stability of boreal forest is shown in Figure 5. The four most influential variables were insects, fire, harvest, and geology. The partial dependence plots are given in Figure 6. The ablation study values in Figure 7; note that any of the RMSE values that are lower than the original model (the “all_boreal” row) perform better when the variable is removed from the model.

Pixel drills for two locations in this study area showing the time series for fPAR and SIWSI, the two-component indices of the forest canopy stability index, and the impact of fire, insect infestations, and logging are shown in Figure 8 and Figure 9.

4. Discussion

The results confirmed that the stability index is tracking both natural and anthropogenic disturbances in the Canadian boreal forests of Quebec and Ontario. Forested areas with low canopy stability index values (i.e., more unstable canopies) were shown to be spatially correlated with recorded logging, fire, and insect disturbances (Figure 5, Figure 6 and Figure 7). The causal interpretation of these disturbances is straightforward: logging removes canopy trees and harvesting in these boreal forests using clear-fell where most canopy trees in the affected area are removed; forest wildfires, above a certain level of severity, combust living canopy foliage scorch tree canopies [60]; and insect infestations can also cause extensive defoliation of tree canopies [61]. Stable forest northwards can be explained in part. First, the growing season and temperatures decrease at higher latitudes, resulting in fewer insect infestations [62]. Second, the logging frontier has moved progressively north, especially since the 1970s and the onset of contemporary forest management practices [3].

Five non-disturbance environmental variables were also shown to be spatially correlated with stability index values, listed in order of importance: geology, TWI, landform, roads, and (non-forestry) land use (Figure 5, Figure 6 and Figure 7). When the human and natural disturbances are removed from the GBM modeling, geology, TWI, and landform are the two more important variables (Figure S8). Previous studies have shown that TWI and landform are related to landscape level gradients in forest flammability and susceptibility to insect infestation, as well as forest site productivity and hence locations more likely to be logged [3,44].

The two examples of the time-series pixel drills illustrate two locations with different disturbance histories in terms of logging, wildfire, or insect infestations. The pixel drill results in Figure 8 had a fire event at the start, which was followed by canopy recovery as indicated by the time series showing the increase in fPAR (greater photosynthesis capacity) and decrease in SISWI (canopy water stress), and an insect infestation towards the end. The pixel drill results in Figure 8 show a stable canopy until the site was impacted by logging. The slope of the CV metrics therefore helps reveal prior instabilities that would not be apparent from a single “snapshot” remotely sensed image.

Our approach provides information that can support decision-makers in understanding the impacts of forest management and natural disturbances on the ecosystem integrity (sensu [63]) of boreal forests. Boreal forests in Canada are currently not experiencing extensive deforestation, however, there is evidence that they are being subject to ecological degradation by current forest management where the primary aim is wood supply [3,64], including how these impacts are now interacting with climate change [2]. While there is no internationally agreed definition of forest degradation, reported impacts include changes in forest structure, reduction in older forests, including old growth, and loss of habitat for disturbance-sensitive species [65]. Forest degradation can be understood as only being the result of the direct impact of human activities such as forest management [63], with wildfires and insect infestation being understood as components of a forest’s natural ecosystem dynamics. However, human-influenced climate change is driving wildfire and insect infestation disturbances outside their natural range, and integrated responses are now required for forest management that factor in how forest management and climate change interact [66]. The time series of the CV slope metrics are useful because they reveal where disturbed forest canopies have recovered, are recovering, or have likely suffered failed regeneration, post-disturbance.

We used 500 m resolution MODIS data to calculate the forest canopy stability because of the 20-year time series, with coverage of the entire Earth’s surface every 1–2 days at resolutions ranging from 250 m to 1 km [67], providing the ability to obtain high-quality, largely cloud-free imagery of fPAR and SIWSI. Comparable results were found by [18] for boreal forests in Siberia, Russia. The approach therefore is likely applicable across the circumpolar boreal forest biome. Ongoing calculation and updating of the stability index is feasible given that continuity plans are in hand to ensure the persistence or functional replacement of each MODIS product beyond the end-of-life of the Terra and Aqua MODIS platforms [68]. In addition, it is also technically possible to generate a comparable time series of the stability index using Landsat and Sentinel-2-sourced data at finer pixel resolutions, though further research is required. However, even as finer resolution stability indices based on time-series analyses become available, the 500 m pixel resolution used here will remain useful for global and regional modeling and conservation assessment purposes. The reason is that ecosystem processes operate across a range of space/time scales [69], and there are impacts on the broader dynamics of the surrounding forest stands and landscape, such as wildlife habitat connectivity [70], that are only evident at larger scales.

A limitation of our analysis was representing the wildfire and harvesting records as simply binary data, i.e., the presence/absence of a wildfire or harvesting event at a pixel. However, these records also contain data on the year disturbances occur and their intensity. By incorporating these additional attributes into the GBM analysis, further insights could be gained into the accumulated impacts of disturbance events on forest stability, and subsequent investigations into the consequences for forest degradation and recovery. A second issue that warrants further research is to investigate the use of other indices in addition to fPAR and SIWSI, including PIS [71] and the Enhanced Vegetation Index (EVI) [72].

While the MODIS data are available to calculate the stability index globally, the extent to which our approach to examining the causal disturbance factors is applicable in other boreal forests depends on the availability of spatial databases on historic harvesting, wildfire, and insect infestations. This is certainly feasible in other Canadian provinces and the USA, as these jurisdictions have comparable systematic forest and disturbance data that are readily available in the public domain [73]. However, this is not the case for other boreal forests [18] and in many developing countries where land use history data are more limited [74].

5. Conclusions

Insights into forest degradation from the stability index require access to ancillary time-series spatial data in order to attribute whether identified disturbed areas are due to human or natural causes. Logging, wildfires, and insect infestations have varying ecological impacts on forest structure, taxonomic composition, and ecosystem processes, including varying rates of post-disturbance regeneration, and insights into causal factors can inform forest management strategies. Furthermore, locations are increasingly experiencing multiple disturbances, and new climate-related pressures with compounding and aggregating impacts on forest ecosystem integrity, requiring more integrated responses. The stability analyses and associated time-series spatial datasets are available via an interactive web-based mapping application, which can be used to help identify ecologically degraded areas in need of restoration and stable forest warranting protection.