Open AccessArticle

Monitoring Changes in the Enhanced Vegetation Index to Inform the Management of Forests

Department of Ecology and Evolutionary Biology, University of Toronto, Toronto, ON M5S 3B2, Canada

Apex Resource Management Solutions, Ottawa, ON K2A 3K2, Canada

Department of Biology, McGill University, Montreal, QC H3A 1B1, Canada

⁴

Quebec Centre for Biodiversity Science, McGill University, Montreal, QC H3A 1B1, Canada

⁵

Group on Earth Observations Biodiversity Observation Network, Montreal, QC H3A 1B1, Canada

Author to whom correspondence should be addressed.

Remote Sens. 2024, 16(16), 2919; https://doi.org/10.3390/rs16162919

Submission received: 25 June 2024 / Revised: 27 July 2024 / Accepted: 6 August 2024 / Published: 9 August 2024

(This article belongs to the Special Issue Remote Sensing Applications for Forest Ecosystem Monitoring and Spatial Modeling)

Download

Browse Figures

Figure 1
The study area is the Algonquin Provincial Park and the surrounding area, which we refer to as the Algonquin Greater Park Ecosystem (AGPE). Elevation (in meters above sea level, masl) is shown in (a). The forest’s mean age in 2019 is shown in (b). Protected areas within the AGPE are shown in (c). About 16% of forest pixels belong to protected areas. The geographic distribution of three disturbance agents in the period 2002–2020 is shown in (d). The gray color represents forested pixels and white non-forest pixels. The dashed black line shows the study area footprint (≈15,000 km2 or 1.5 M ha). The solid black line shows Algonquin Provincial Park’s footprint. The two perpendicular dashed lines in (a) divide the study area into quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the description of the results. Map projection is NAD83 Statistics Canada Lambert, EPSG: 3347. "> Figure 2
Main methodological steps used in this study. After downloading MODIS EVI data, pixels were filtered to only keep pixels with good quality data, within forested areas. EVI time series were then created and processed with three BFAST algorithms: bfast, bfast01, and bfastclassify. The bfast algorithm decomposes a time series into a seasonal component, trend, and noise. Using piece-wise linear regression on the trend component, it detects one or more breaks (if any). Here, we use three main outputs provided by bfast: type of break, magnitude of break, and time of break with 95% confidence intervals (CIs) The bfast01 algorithm runs a seasonally adjusted regression model on the ts and only detects the major break (if any). The bfastclassify algorithm then uses bfast01’s output to classify trends into one of eight possible trend types (<a href="#remotesensing-16-02919-f0A1" class="html-fig">Figure A1</a>). Only the magnitude of break values estimated by bfast was used in boosted regression trees (XGBoost models) to explore their relationship with predictor variables (dashed box at bottom, Equation (<a href="#FD1-remotesensing-16-02919" class="html-disp-formula">1</a>) in main text). Satellite icon from <a href="http://flaticon.com" target="_blank">flaticon.com</a> (accessed on 29 April 2022). "> Figure 3
Spatial distribution of EVI negative and positive breaks in the AGPE from 2003 to 2022. Most of the breaks were found in the eastern half of the AGPE and more particularly in the northeast quadrant. There were 11,871 pixels with negative breaks (red pixels) and 3893 pixels with positive breaks (cyan pixels). These breaks were estimated with the bfast algorithm. The dashed black line shows the study area footprint. The solid black line shows Algonquin Provincial Park’s footprint. The two perpendicular dashed lines divide the study area into four quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the interpretation of results. Map projection is NAD83 Statistics Canada Lambert, EPSG: 3347. "> Figure 4
Spatial distribution of EVI trend types in the AGPE from 2003 to 2022. These trends were produced with the bfastclassify algorithm. The trend types are those proposed by de Jong et al. [<a href="#B57-remotesensing-16-02919" class="html-bibr">57</a>]. Abbreviations—MIG: monotonic increasing, greening trend (n = 33,683); MDB: monotonic decreasing, browning trend (n = 4981); IInb: interruption, increasing trend with a negative break (n = 11,637); RBG: reversal, browning to greening trend (n = 11,654). (Trends MIGpb, MDBnb, IDpb, and RGB are not shown given their low percentages). The dashed black line shows the study area footprint. The solid black line shows Algonquin Provincial Park’s footprint. The two perpendicular dashed lines divide the study area into four quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the interpretation of results. Map projection across all figures is NAD83 Statistics Canada Lambert, EPSG: 3347. "> Figure 5
Predictors of the magnitude of EVI breaks (negative breaks in red and positive breaks in cyan tone boxes) from 2003 to 2022. The ranking reflects feature importance using the gain metric as estimated by XGBoost models. The same five predictors are in the top five but with slightly different rankings (connecting lines with slopes) except for the summer climate moisture index with a 3-year lag, which ranks fourth in both (connecting line with no slope). Forest protection status is low-ranking for both types of breaks. The XGBoost models were run with subsets of the detected breaks (negative breaks, n = 116 records; positive breaks, n = 3263 records). "> Figure A1
Schematic of break and trend types employed in this study. A time series (ts) can be characterized by the presence or absence of breaks and trends. Breaks represent abrupt changes in a ts whereas trends represent gradual changes. Here, we refer to three major groups of trends: monotonic (a,b), interruption (c), and reversal trends (d). Similarly, we refer to two break types: negative (red downward arrow) and positive (green upward arrow) breaks. Trends can be monotonic, either increasing (green slopes) or decreasing (orange slopes). The former are referred to as greening trends and the latter as browning trends. Monotonic trends can show no breaks (a) or show one or more breaks (negative or positive) but in concordance with the slope of the trend segments (b). Conversely, interruption (c) and reversal (d) trends are characterized by having a break type in discordance with the slope of the trend segments. Interruption trends can have two positive trend segments divided by a negative break and vice versa (two negative trend segments divided by a positive break). Reversal trends have opposite trend segments divided by a negative or positive break. Lastly, some ts may not change or show changes that are too small to be detected with the methods employed (horizontal gray dashed line in (a)). Here, we use the trend classification proposed by de Jong et al. [<a href="#B57-remotesensing-16-02919" class="html-bibr">57</a>]—MIG: monotonic increasing, greening trend (without breaks) (bottom line in (a)); MDB: monotonic decreasing, browning trend (without breaks) (top line in (a)); MIGpb: monotonic increasing, greening trend with a positive break (top set of lines in (b)); MDBnb: monotonic decreasing, browning trend with a negative break (bottom set of lines in (b)); IInb: interruption, the increasing trend with a negative break (top set of lines in (c)); IDpb: interruption, decreasing trend with a positive break (bottom set of lines in (c)); RGB: reversal, greening to browning trend (top two sets of lines in (d)); RBG: reversal, browning to greening trend (bottom two sets of lines in (d)). "> Figure A2
Density plots of forest EVI magnitude of breaks in the AGPE from 2003 to 2022. The number of breaks and their magnitudes broken down by year are shown. Extreme magnitude values have been omitted to aid visualization. Vertical lines show medians—solid: yearly; dashed: entire time series. The total number of breaks was 15,904 (11,969 negative and 3935 positive). The time of break was rounded up prior to plotting which caused 2003 breaks (8 negative and 18 positive) to be part of 2004. No breaks were detected in 2022. "> Figure A3
Ranking of all predictors (features) used in the XGBoost models. Panel (a) shows predictors of negative break magnitudes and panel (b) shows those of positive break magnitudes. These models were run with subsets of the detected breaks (negative breaks, n = 116 records; positive, n = 3263 records). Variable abbreviations—for_age: forest age; dd5_wt: winter degree days above 5 °C; cmi_sm: summer climate moisture index; for_con: percentage of conifers; for_pro_0: non-protected forest; lag#: 1-, 2- or 3-year lags. The protected forest variable is not present in (a) given its lack of importance in explaining negative breaks. "> Figure A4
Partial dependence plots of magnitude of negative break predictors (features). The relationships between predictors and response (magnitude of EVI breaks) variables were mostly non-linear. Plots were created from the output of XGBoost. The model was run with a subset of all detected negative breaks (n = 116). The values on the y-axes are absolute values of negative magnitudes. Variable abbreviations—for_age: forest age; dd5_wt: winter degree days above 5 °C; cmi_sm: summer climate moisture index; for_con: percentage of conifers; for_pro_0: non-protected forest equals 1, protected equals 0; lag#: 1-, 2- or 3-year lags. "> Figure A5
Partial dependence plots of magnitude of positive break predictors (features). The relationships between predictors and response (magnitude of EVI breaks) variables were mostly non-linear. Plots were created from the output of XGBoost. The model was run with a subset of all the detected positive breaks (n = 3263). Variable abbreviations—for_age: forest age; dd5_wt: winter degree days above 5 °C; cmi_sm: summer climate moisture index; for_con: percentage of conifers; for_pro_0: non-protected forest equals 1, protected equals 0; lag#: 1-, 2- or 3-year lags. "> Figure A6
Geographic distribution of trends in Algonquin Park and the surrounding area, which we refer to as the Algonquin Greater Park Ecosystem (AGPE). All maps show trends that were derived from the output of bfastclassify. Compared to greening trends (MIG) which occurred throughout the AGPE (a), browning trends (MDB) mostly occurred in the NE quadrant (b). Most increasing trends with negative breaks (interruptions, IInb) occurred in the NW quadrant (c) while most of the relatively few decreasing trends with positive breaks (interruptions, IDpb) occurred in the NE quadrant (c). Notably, browning to greening reverse trends (RBG) co-occurred with browning trends in the NE quadrant (d). The dashed black line shows the study area footprint (1.5 M ha). The solid black line shows Algonquin Provincial Park’s footprint. The perpendicular dashed lines divide the study area into four quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the description of the results. Map projection across all figures is NAD83 Statistics Canada Lambert, EPSG: 3347. ">

Versions Notes

Abstract

In the absence of forest ecosystem time series data, monitoring proxies such as the enhanced vegetation index (EVI) can inform the capacity of forests to provide ecosystem services. We used MODIS-derived EVI at 250 m and 16-day resolution and Breaks for Additive and Seasonal Trend (BFAST) algorithms to monitor forest EVI changes (breaks and trends) in and around the Algonquin Provincial Park (Ontario, Canada) from 2003 to 2022. We found that relatively little change occurred in forest EVI pixels and that most of the change occurred in non-protected forest areas. Only 5.3% (12,348) of forest pixels experienced one or more EVI breaks and 27.8% showed detectable EVI trends. Most breaks were negative (11,969, 75.3%; positive breaks: 3935, 24.7%) with a median magnitude of change of −755.5 (median positive magnitude: 722.6). A peak of negative breaks (2487, 21%) occurred in the year 2013 while no clear peak was seen among positive breaks. Most breaks (negative and positive) and trends occurred in the eastern region of the study area. Boosted regression trees revealed that the most important predictors of the magnitude of change were forest age, summer droughts, and warm winters. These were among the most important variables that explained the magnitude of negative (R² = 0.639) and positive breaks (R² = 0.352). Forest composition and protection status were only marginally important. Future work should focus on assessing spatial clusters of EVI breaks and trends to understand local drivers of forest vegetation health and their potential relation to forest ecosystem services.

Keywords:

abrupt forest change; gradual forest change; forest disturbance; forest ecosystem services; BFAST

1. Introduction

Human activities induce changes in forest disturbance regimes that can convert once vast carbon sinks, like those of Canada’s forests, into carbon sources [1,2]. In Canada, an average of 27 Mt of C year⁻¹ was released into the atmosphere due to fires (corresponding to the burning of 2 M ha of forests) from 1959 to 1999 [3]. Then, in 2023, the Canadian forest fire season was record-breaking, burning more than 18 M ha of forests (about 6.7 times the prior 10-year average). Similarly, insect disturbance affects vast areas of Canadian forests not only reducing forest productivity but also making forests more prone to fire [4]. Furthermore, droughts not only kill trees outright but also cause tree stress, which impacts productivity negatively, reducing forest ecosystem services [5,6,7]. Importantly, droughts can make trees more susceptible to insect damage [8] and show multi-year temporal lag effects [6,9]. All these forest disturbances are (and will continue to be) exacerbated by climate change [10,11], with important impacts on the forests of northern latitudes [12,13,14].

Given the changes in the disturbance regimes, forest ecosystem services need to be “monitored for action” as is being advocated by the biodiversity monitoring community [15]. That is, the information gained through monitoring needs to be translated into actions at the level of policy and on the ground [16]. There is no single monitoring approach that can provide all the information needed for effective action, but approaches that adopt principles of adaptive monitoring [17] and change detection-attribution frameworks [18,19] are likely to be more effective. Moreover, monitoring approaches that leverage open data are more likely to succeed and scale up globally [20]. Furthermore, it is increasingly important to monitor not only the ecological supply of forest ecosystem services but also other dimensions of ecosystem services [21,22] as well as the drivers affecting these services.

In the absence of accessible forest aboveground biomass high-resolution time series, remote sensing-derived vegetation indices (VIs) can be useful proxies of forest carbon services and other forest ecosystem services. Commonly used VIs are readily and openly available online and have become important variables for ecological analyses [23,24,25]. These VIs, which are available at relatively fine spatiotemporal resolutions for periods of a decade or more, are good indicators of forest health status as they reflect the effects of environmental and ecological factors [26,27]. Indices such as NDVI (normalized difference vegetation index) and EVI (enhanced vegetation index), which measure vegetation greenness, are common inputs for the estimation or validation of forest greenness, biomass, and carbon [28,29,30,31]. Decision-makers are often in need of simple and easy-to-interpret values or thresholds of ecological change to reach consensus and prompt action across government levels. Easy-to-implement monitoring approaches that leverage easy-to-interpret VIs can help translate knowledge to action, which is key for achieving biodiversity conservation goals [15,16] as well as managing ecosystem services.

The analysis of EVI time series of forested ecosystems can provide important insights about EVI dynamics and forest health. Negative abrupt changes in the EVI time series (known as trend breakpoints or simply breaks) can be due to pulse disturbances (e.g., intense forest fires or windthrows [32]) while positive abrupt changes can be due to boosts of natural regeneration or reforestation (e.g., plantation). Then, gradual change (i.e., trends) in a time series can be due to press disturbances [33] like extended droughts (negative trend) or sustained good growing conditions (positive trend). Thus, breaks and trends can be either positive or negative (Figure A1) depending on the disturbing or enhancing factor at play. The quantification of breaks and trends in the forest EVI time series and their link with climatic and environmental drivers can help, therefore, to determine the current and future state of forested landscapes. A forest that is greening as opposed to one that is browning is likely healthier and hence, provides more ecosystem services.

To quantify EVI changes and assess forest health in Algonquin Provincial Park and the surrounding area, we use remote sensing data and change detection algorithms. Specifically, we focus on the following objectives: (1) detect and quantify significant breaks in the forest EVI, (2) detect and quantify significant trends in the forest EVI, and (3) explain variations in the magnitudes of forest EVI breaks through the forest- and climate-related predictors. We predict that (1) protected forests will show more positive breaks and greening trends (but less browning) compared to non-protected forests as they experience less anthropogenic pressures; (2) negative abrupt changes in the EVI are more related to disturbances and climate stress compared to positive changes, which are likely more related to favorable climate conditions; and, (3) forest age and composition are important for understanding EVI abrupt change as drivers affect trees differently depending on these factors [34,35].

2. Materials and Methods

2.1. Study Area

The study area encompasses the Algonquin Provincial Park (7630 km², 763,000 ha) and the surrounding area (hereafter, Algonquin Park and Greater Park Ecosystem (AGPE)), covering a total of about 15,000 km² or 1.5 M ha (Figure 1). The area is located in south-central Ontario (Canada), on Precambrian Canadian Shield bedrock, and contains about 15% water and 2500 lakes [36]. Elevation varies from about 150 m above sea level (masl) in the east to about 600 masl in the west (Figure 1a) [37]. This gentle altitudinal gradient has important climatic consequences—the west is relatively cooler and wetter than the east [37]. Most of the older forests are located in the eastern part of the AGPE (Figure 1b) and about 16% of the area is protected (Figure 1c) with little or no timber harvesting activities.

The AGPE has seen its share of disturbances. Most of Algonquin Park has burned at least once over the past 150 years [36] and timber harvest and insect pest outbreaks have occurred periodically in the area (Figure 1d). Active forest and fire management (fire suppression and prescribed burning) take place in the study area [36]. The AGPE can be considered a working landscape [38] that has been impacted by natural disturbances and human activities. This makes the AGPE a great case study to track forest vegetation change across time and space and shed light on the implications of this change for forest ecosystem services.

2.2. Data

2.2.1. Vegetation Index

We used the Moderate Resolution Imaging Spectroradiometer (MODIS) enhanced vegetation index (EVI) at 16-day and 250-m resolution for the period 1 January 2003 to 31 December 2022 (20 years) as the main response variable (Table A1). EVI data (product version: MOD13Q1 v061, Terra satellite) were downloaded from the Land Processes Distributed Active Archive Center’s (LP DAAC) server using the R statistical programming language (version 4.1.1, [39]) and the MODIStsp package (version 2.0.9, [40]). We filtered out poor-quality pixels based on the approach suggested by Samanta et al. [41]. Furthermore, non-forest pixels were masked out with a forest land cover data layer for the year 2003 by Hermosilla et al. [42]. MODIS EVI values have a valid range from −2000 to 10,000 while forested landscapes are mostly associated with values > 1800.

2.2.2. Forest Attributes and Disturbance

We used forest age [43,44], forest composition [43,45], and forest land cover [42] (Table A1). Details on the development and validation of these datasets have been described elsewhere [42,44,45,46]. With respect to disturbance data, we used forest fires and harvest [46] as well as Canada Landsat Disturbance (CanLaD, [47]) and Hansen et al.’s [48] Global Forest Change (GFC) 2000–2022 product (version 1.10). Furthermore, we used the forest insect damage event database (polygon-based) provided by the Government of Ontario [49].

2.2.3. Climate and Other Data Sources

We extracted 30 climate-related variables from ClimateNA (version 7.31, [50]) for the period 1998–2021 and converted them to gridded data (i.e., raster formatted data). We linked these gridded data to our response data (i.e., EVI breaks, see Section 2.3.1 below) by matching pixel ID and year. Then, we created time-lagged versions of these variables by keeping the pixel ID match but modifying the year match by 1, 2, and 3 years (1-, 2-, and 3-year temporal lags). We applied variance inflation factors (VIF, threshold: VIF < 5) to the full set of climate variables, keeping only six climate variables that were not collinear: winter degree days above 5 °C (winter DD5) with a 1-year lag, winter DD5 with a 3-year lag, summer climate moisture index (summer CMI), summer CMI with a 1-year lag, summer CMI with a 2-year lag, and summer CMI with a 3-year lag. Warm winters and hot dry summers can reduce water availability and increase evapotranspiration, which in turn affect tree growth, health, and the provisioning of forest ecosystem services [5,51,52].

We also used Canada’s protected areas boundaries [53] and elevation data (NASA Shuttle Radar Topography Mission Global, 1 arc second ≈ 30 m resolution, product: SRTMGL1 v003). A list of key data sources is shown in Table A1 and more details are provided in Appendix A.1.

2.3. Statistical Methods

2.3.1. EVI Change Detection

We used Breaks for Additive Seasonal and Trend (BFAST, [54,55]) statistical methods to detect abrupt (i.e., breakpoints or breaks) and gradual changes (i.e., trends) in the forest EVI across time series pixels while accounting for seasonal variation. (Figure A1 in Appendix A.2 shows a schematic of the different types of breaks and trends mentioned here). BFAST’s algorithms have been used in several studies analyzing vegetation dynamics [56,57,58,59,60,61]. We used the BFAST R package (version 1.6.1) and three algorithms-standard bfast (hereafter bfast), bfast01 and bfastclassify.

Thus, bfast decomposes a time series into seasonal, trend, and remainder components and integrates these with methods for detecting and characterizing time series changes [54,55,62]. Trend breakpoints are detected using piece-wise linear regression over segments of the detrended time series. The bfast function finds one or more trend breakpoints (hereafter, breaks) in a pixel time series. Among the key outputs, bfast provides the magnitude, direction, or type (positive or negative) and the approximate timing of each break (together with 95% confidence intervals, hereafter 95% CIs (Figure 2).

Furthermore, we ran the bfast01 algorithm, which provides input to the bfastclassify algorithm to categorize trends (Figure 2). Contrary to the standard bfast algorithm, bfast01 runs a seasonally adjusted regression model to find the single biggest break (if any) across the entire time series. If a statistically significant break is found, then the algorithm fits two trend segments based on piece-wise linear regression, one on each side of the break. If no statistically significant break is found, the algorithm fits a single trend line to the time series. Each pixel time series processed by bfast01 is classified by bfastclassify into a trend type as proposed by de Jong et al. [57].

One of the key parameters for BFAST’s algorithms is h, which determines the temporal window size between potentially detected breaks. h should be set according to the frequency and magnitude of the particular disturbance events and vegetation responses attempted to be captured [60,63]. Often, these details are not known in advance. Here, we set h to 0.05, which is equivalent to one year of data given our time series (23 time steps per year when using the MODIS EVI 16-day temporal resolution). This choice of h means that once a break is detected, a new break can only be detected after one year of observations has passed. A sensitivity analysis showed that h = 0.05 provided the best compromise to leverage the disturbance data used here and produce annual summaries (Table A2 and Table A3).

After prepossessing and filtering out poor-quality pixels (e.g., pixels with high cloud cover, high aerosols, possible shadow, etc.), we ran BFAST algorithms on two decades of forest EVI time series (2003–2022). Analyses were run using a 0.01 significance level for bfast and bfast01 and 0.05 for bfastclassify. All BFAST analyses presented were conducted using h = 0.05. (For sensitivity analyses of h, see Table A2 and Table A3). To attribute forest EVI negative breaks to disturbances, we overlapped detected break and disturbance pixels in space and time. We assume that an overlap meant that a particular disturbance was linked to a particular EVI break.

2.3.2. Modeling EVI Magnitude of Change

To determine the relative contribution of climate and forest-related variables in explaining EVI magnitude of change we used boosted regression trees (BRT [64]). BRT is a statistical (machine) learning technique that is better suited than traditional statistical approaches to deal with non-linearities and interactions that are likely present in ecological data [64]. We used the Extreme Gradient Boosting (XGBoost) library to run BRT in an efficient, flexible, and scalable manner [65].

The following model formulation was used with XGBoost:

d Y = f (X 1, X 2)

(1)

where dY is a vector containing the magnitude and sign (positive or negative) of EVI breaks. X1 represents a matrix of six climate-related predictors: winter DD5 with a 1-year lag, winter DD5 with a 3-year lag, summer CMI, summer CMI with a 1-year lag, summer CMI with a 2-year lag, and summer CMI with a 3-year lag. X2 represents a matrix of three forest-related predictors: forest age, percent conifer trees, and forest protection status. Each row in the data frame corresponds to the magnitude of a distinct break (one or more breaks per time series pixel are possible) linked in time and space with predictor variables (based on pixel ID and year).

A total of nine predictors were used in each model (Table 1). We ran one model for negative breaks and another for positive ones. XGBoost (version 1.7.3) was computed using the Python programming language (version 3.10). More details about the steps taken are given in Appendix A.3. Figure 2 summarizes the major steps described above with a high-level methodological flow chart.

3. Results

3.1. EVI Changes in Time

The study area contained 241,378 forest pixels at 250 m resolution. Around 4% (n = 10,472) of these could not be processed by BFAST’s algorithms mainly due to data gaps, leaving 230,806 pixels (12,416.83 km², 1,241,683 ha). A total of 15,904 breaks (11,969 negative and 3935 positive, 75.3% and 24.7%, respectively) were found by the bfast algorithm (Table 2). Figure A2 shows the temporal distribution of negative and positive EVI breaks. A peak of negative breaks (n = 2487; 21% of all negative breaks) occurred in the year 2013 while no clear peak was seen among positive breaks over the two decades analyzed. The median magnitude of change for negative breaks was −755 while that for positive breaks was 723. The strongest median negative break (−1141) occurred in 2004 and the highest median positive break (1277) occurred in 2018. Comparing the 20-year median with the annual median, five consecutive years (2018 to 2022) were more negative than the overall negative median. Among positive breaks, four consecutive years (2008 to 2011) showed consecutive medians that were lower than the 20-year positive median.

3.2. EVI Changes Throughout the Study Area

Most forest pixels that showed breaks only had one break (74.9%) over the two decades while some showed more than one break. Of the latter pixels (>1 break, 25.1%), most had only one type of break (negative or positive) while a few (3.3%) showed both negative and positive breaks over the two decades analyzed. After collapsing multiple temporal breaks by pixel, we found that 5.3% (n = 12,348) of all forest pixels experienced EVI breaks (abrupt changes) over the 20-year study period (Table 2). On average, in any given year there was a 0.34% chance that a forest pixel would experience a break in the study area. Most forest EVI breaks (both negative and positive pixels) occurred in the eastern region of the AGPE (Figure 3). There were more negative break clusters in the northeast quadrant compared to other quadrants.

Comparing breaks by forest protection status, 95.4% of all breaks occurred across non-protected areas, which represent 84.1% of the study area, and 4.6% of breaks occurred across protected areas, which represent 15.9% of the study area. Negative break pixels amount to 4.8% of all non-protected forest pixels and only 1.0% of all protected forest pixels (Table 2). Positive breaks amount to 1.5% of all non-protected forest pixels and 0.5% of all protected ones. Positive breaks were bigger in protected areas (median = 794.5, mean = 819.6, n = 256) compared to non-protected forests (median = 719.4, mean = 801.5, n = 3679, Wilcoxon test = 510,702, p-value = 0.0236). With respect to negative breaks, the difference between protected (median = −765.6, mean = −854.3, n = 470) and non-protected forests (median = −754.6, mean = −829.5, n = 11,499) was not significant.

Attribution of negative breaks to disturbance data (remote sensing derived) was conducted with varying degrees of success depending on the data source used (Table A3). For instance, only 16% of detected negative breaks could be matched (in space and in time using time of break 95% CIs) with forest harvest pixels, and less than 1% of breaks could be matched with insect damage (moderate and severe tree infestation and mortality) (h = 0.05 heading in Table A3). When using CanLaD (fire and harvest), a 28% match with negative breaks was achieved and with GFC (a disturbance dataset with no labels) the matching went up to 60%. Thus, in the best of cases, about 40% of the negative breaks found in this study could not be matched to disturbance data. This suggests that many of the negative EVI breaks are likely associated with low or moderate non-stand replacing disturbances or other factors (i.e., soil conditions) that left a strong enough signal that was detected with EVI.

The bfastclassify algorithm revealed that 27.8% of forest pixels had statistically significant trends (Table 3). Similarly, about 14.6% of forest pixels experienced greening (monotonic increasing trends without breaks, MIG) while 2.2% experienced browning (monotonic decreasing trends without breaks, MDB). About 5% of forest pixels showed interruption trends (increasing trend with a negative break, IInb) and 5% showed reversal trends (browning to greening, RBG). About 2.6% of forest pixels showed interruption trends characterized by a decreasing trend with a positive break (IDpb). Breaking down the results by protection status, 15% of all non-protected forest pixels showed greening trends (MIG) compared to 8.4% of protected forest pixels. Similarly, MDB trends represent 2.2% of trends in non-protected areas compared to 1.6% in protected ones. Furthermore, 5.3% of all non-protected forest pixels showed IInb trends while only 2.4% of protected forest pixels showed this trend. Likewise, 5.5% of all non-protected forest pixels showed RBG trends compared to 1.4% of all protected forest pixels. That is, non-protected forests showed about twice the percentage of IInb trends and about four times the percentage of RBG trends compared to protected forests (5.3% vs. 2.4% and 5.5% vs. 1.4% respectively). All other trend types including monotonic increasing trends (greening) with positive breaks (MIGpb), monotonic decreasing trends (browning) with negative breaks (MDBnb), and interruption trends characterized by decreasing trends with positive breaks (IDpb) represented <1% of all forest pixels (Table 3).

In terms of spatial distribution, similar to the breaks, most of the statistically significant trends occurred in the eastern half of the AGPE (Figure 4). Browning trends were mostly detected in the northeast quadrant while greening trends were more evenly distributed. There was a noticeable concentration of interruption trends (type IInb) in the northwest quadrant of the AGPE.

3.3. Drivers of Change

Most detected breaks had a wide time of break 95% CI range (negative: min = 1 month, max = 82 months or 6.8 years, average = 11 months; positive breaks: min = 1 month, max = 74 months or 6.2 years, average = 8 months). Hence, we subset the data to keep records that were more likely to show a stronger signal of abrupt change. We ran XGBoost on a subset of negative breaks with narrow ranges (<1.5 months) and with a magnitude of change > |−500| (n = 116). For positive breaks, we focused on all breaks with a magnitude of change > 500 (n = 3263), regardless of the time of break 95% CI range. For positive breaks, we also removed breaks that overlapped with GFC disturbance data. The improvement of our fine-tuning was reflected in the larger R² adjusted values for the models that used the data subsets compared to those that used all the records. (Table A4 in Appendix A.3 provides the results for the different models that were run).

XGBoost models explained 64% of the variance among negative breaks (R²-adj = 0.639, RMSE = 292.5) and 35% among positive breaks (R²-adj = 0.352, RMSE = 310.4). Forest age, winter DD5 with a 1-year lag, and summer CMI with 1 to 3-year lags were among the five most important variables for negative and positive breaks but with slightly different rankings (Figure 5). For negative breaks, forest age was the most important predictor while it was the second most important for positive breaks. Winter DD5 with a 1-year lag was the most important predictor for positive breaks and the second most important for negative breaks. The effect of forest age was about twice that of winter DD5 (1-year lag) among negative breaks but about the same among positive ones (Figure A3a). Summer CMI ranking was slightly different, whereas, for negative breaks, 1-year lag summer CMI ranked third, it ranked fifth for positive breaks, and vice versa. Summer CMI with a 3-year lag ranked fourth for both types of breaks. While the percentage of conifer composition ranked seventh among negative breaks, it was last (ninth) among positive ones. The predictive strength of forest protective status was relatively very small for positive breaks (ranking second to last) and not at all important for negative breaks (being absent in Figure A3a).

BRT partial dependence plots for the magnitude of breaks showed interesting patterns (Figure A4 and Figure A5 Appendix A.3). The magnitude of change of negative breaks was high (∼2100) and remained so for younger forests until about 90 years when it dropped sharply (Figure A4a). The percentage of conifer composition showed a similar curve for negative breaks, starting at a higher EVI value (∼2300) and showing a more gradual decrease beginning at around 40% (Figure A5g). In contrast, positive breaks showed decreasing non-monotonic curves (starting at ∼850) with forest age but with two small bumps, one at about 60 years and a second one at about 90 years of age (Figure A5b). With respect to the percentage of conifers, positive breaks showed a non-monotonic curve with no clear pattern (Figure A5i).

Winter DD5 with a 1-year lag showed a high negative change curve (∼2000) throughout its range (Figure A4b). For positive breaks, this predictor increased non-linearly with warmer winters and peaked at 850 around 7-degree days (Figure A5a). Furthermore, the 3-year lag version of this predictor showed a similar pattern for negative breaks with a constant high value around 2100 (Figure A4h) but for positive breaks, it showed a relatively flat non-monotonic curve that fluctuated around 800 (Figure A5f). This means that the lagged effects of warmer winters on abrupt brownness (negative breaks) were more than twice its effect on abrupt greenness (positive breaks).

Summer CMI predictors had similar curves for negative breaks; that is, they all showed high constant negative change values (between 2000 and 2200) throughout much of their range with sharp drops at one or both extremes (Figure A4c–f). Interestingly, for positive breaks, summer CMI predictors showed lower non-monotonic curves (Figure A5c–e,g). Notably, the curve of summer CMI with a 2-year lag was roughly U-shaped. It showed relatively high values (above 950) associated with drought-like conditions, corresponding to negative CMI values. There was a dip when CMI values were close to zero and an increasing trend for values higher than 5 (Figure A5c).

When using the data subsets, the effect of protection status on the magnitude of positive abrupt changes (positive breaks) is relatively minor, as shown by the shallow increasing slope between protected and non-protected forests (Figure A5h). Interestingly, non-protected status was associated with a slightly higher value of greenness (positive break magnitude) compared to protected status. This may be due to plantation regrowth occurring in timber harvest areas. Protection status has no association with the magnitude of negative abrupt changes (Figure A4i).

4. Discussion

Our goal was to gain a better understanding of the relationship between EVI breaks (abrupt changes), trends (gradual changes), and forest health. We shed light on this relationship by quantifying forest EVI breaks and trends in the AGPE and explaining variation in break magnitude with forest and climate-related predictors.

Forest vegetation changes arise in part from stand-replacing (e.g., forest fires) and non-stand-replacing disturbances (e.g., hail tree damage). Negative EVI breaks likely reflect the impact of stand-replacing disturbances whereas negative trends (browning trends) are likely associated with non-stand-replacing disturbances. Conversely, positive breaks may reflect growth spurts associated with improved climate conditions, competitive release, or successful reforestation activities. Positive trends (greening trends) may reflect sustained growth or recovery after a disturbance. These important forest dynamics can be hidden when assuming that change occurs monotonically; hence, detecting breaks and differentiating between short- and long-term changes allow for a robust way of analyzing forest change [54,56,57,63].

4.1. Changes in EVI

Of the forest pixels having breaks, most of them only had one break (74.9%). This means that most negative breaks were not followed by a detectable positive break or a second negative break. Thus, it seems that recovery from negative breaks is primarily gradual with few post-disturbance bursts of growth. Alternatively, the lack of detection of multiple breaks may be due to the h parameter value used. A smaller value of h (0.025 or six months) detected more breaks per pixel but much fewer breaks (Table A2).

Most of the forest EVI pixels that showed a change in the AGPE were located in non-protected areas (95.6% of all breaks and 91.1% of all trends, respectively). This is in line with the findings by Sharma et al. [66] who found that protected forests change less compared to nearby non-protected forests mainly due to the lack of human activities taking place therein. Also, we found that there are proportionally more positive breaks in protected forest pixels compared to non-protected ones (35.3% vs. 24.2%, respectively), which implies that on average protected forest pixels show more multiple growth bursts per positive break pixel compared to non-protected forest pixels. Contrary to our predictions, there were fewer greening trends among protected forest pixels (but fewer browning trends) compared to non-protected forest pixels, occurring at half the rate of the latter (8% vs. 15.0%). This higher percentage of greening in non-protected forests is likely the result of regrowth from harvesting events or other anthropogenic activities that took place prior to 2003. On the other hand, the low percentage of greening pixels among protected forest pixels may be the result of the low levels of change experienced by these forests, which may be too low to be detected with the approach used here. Alternatively, the lower number of greening trends seen may be due to the EVI saturation effects of denser canopies [67] that are likely found in protected forests. Overall, given that natural disturbances are expected to occur in all forests regardless of protection status [68] and that non-protected forests experience higher rates of anthropogenic disturbance, protected forests in the AGPE are likely better at providing ecosystem services compared to non-protected forests.

In terms of the geographical pattern, the northeast quadrant showed the greatest vegetation activity compared to other quadrants (both in terms of breaks and trends) (Figure 3 and Figure 4). The northeast quadrant also showed clustering of pixels showing long-term browning, greening, and reverse (browning to greening) trends (Figure A6a,b in Appendix A.4). Such results may suggest that the northeast quadrant is recovering from a drought event.

In line with our initial predictions, the magnitudes of positive breaks were bigger in protected forests compared to non-protected ones. Furthermore, in terms of the overall proportion of positive breaks, protected forest pixels show more positive breaks per pixel than non-protected forest pixels. This could imply that on average protected forests show stronger and multiple growth bursts compared to non-protected forests. Interestingly, the difference in the magnitude of negative breaks between protected and non-protected forests was not significant. This is surprising given that most negative breaks and stand-replacing disturbance pixels were found in non-protected forests. However, this could mean that most negative breaks are the result of natural disturbances (e.g., droughts), which are likely to affect all forests within the AGPE.

4.2. Attribution of Change

Once EVI breaks have been detected, linking these breaks to specific drivers becomes crucial but challenging. Stand-replacing disturbances such as intense forest fires and clear-cut harvesting are usually easier to identify with remote sensing data due to their abrupt change signature [33]. Non-stand replacing disturbances such as tree thinning, low to moderate insect or pathogen damage, and low to moderate drought are likely to show more gradual changes and are relatively harder to detect [33]. The former type of disturbance leaves a noticeable change (e.g., negative EVI breaks of large magnitude) in the forest while the latter may leave a more subtle change (e.g., negative EVI breaks of smaller magnitude). Furthermore, some forest breaks may not be caused by a single intense short-lived disturbance event that occurred abruptly (pulse event, [32]) but rather from the cumulative effects of several interacting disturbance factors of varying intensities that have occurred over extended periods of time.

Although timber harvest is an ongoing activity in the non-protected forests of the AGPE, only 16% of abrupt negative changes matched harvest pixels (unmatched harvest pixels 82%). Similarly, only about 1.2% of break pixels matched forest fires (unmatched fire pixels 77%). The low match between negative breaks and harvest pixels could be related to the type of harvesting taking place. Clearcut systems are likely to produce more forest EVI breaks than selection and shelterwood systems, which are commonly used in the study area [37]. Perhaps, the latter two harvesting systems are harder to detect using MODIS EVI data at 250 m resolution, as done here. Very few fire events are likely to occur in the AGPE due to fire management activities. Yet, the relatively low number of EVI break pixels overlapping fire pixels was not expected given that fire disturbance usually produces a distinct abrupt change on the ground.

Less than 1% of EVI negative breaks matched insect damage (only 7.7% of insect damage polygons were intersected). This little overlap with insect damage data is harder to explain. Sommerfeld et al. [9] found that insects can produce high-severity disturbance in temperate forests as noted by on-the-ground local experts. Here, there were very few matches between negative breaks and insect damage. It is possible that the severity of insect damage is diluted at the pixel resolution used (250 m), producing more gradual changes rather than abrupt ones. As such, these gradual changes may not be detected as trend breaks but rather as browning or reversal trends (greening to browning).

When using CanLaD (which captures mostly fire and harvest disturbance pixels) and GFC disturbance data (not labeled), about 28.1% and 60.4% of EVI negative breaks overlapped with these disturbance pixels. Notably, only 13.7% of CanLaD and 10.3% of GFC disturbance pixels overlapped with EVI negative breaks. This means that most disturbance pixels in these two data sources were not linked to negative breaks. (A match among different disturbance datasets (Table A6) in Appendix A.5 shows important mismatches too).

The relatively poor match between EVI negative breaks and the different disturbance data used may reflect the fact that EVI is a one-dimensional metric that may not capture the full complexity of disturbance patterns. Yet, our approach likely detected low-intensity disturbances (e.g., hail, short-lived climate stress) and/or forest degradation [69].

XGBoost models revealed predictors causing variation in the magnitude of breaks as well as strong non-linearities and interactions among predictors. Forest age, winter DD5, and summer CMI were among the top five predictors for negative and positive breaks but with different rankings. Notably, four out of the top five climate predictors ranged in temporal lags from 1 to 3 years. These results are in line with other forest studies that have shown that 2- and 3-year climate lag effects characterize low and moderate disturbance activity clusters in temperate forests [9].

Droughts can induce forest stress and even mortality at large scales that in turn have important effects on forest carbon and global climate [5,7,52,70,71,72]. Our findings hint at the importance of summer droughts (negative CMI values) and warm winters (high winter DD5 values) and their temporal lag effects for understanding the variation in the magnitude of EVI breaks in the AGPE. Our results show that the effects of warm winters and summer droughts can start several years before (3-year lagged summer CMI and 1-year lagged winter DD5) and build up over consecutive years (3-year lag to no-lag summer CMI) prior to the detection of forest EVI breaks. Heat and drought early in the growing season or over consecutive years can have negative effects on forest productivity [35]. Given that global warming is expected to exacerbate the frequency and intensity of droughts around the globe with important consequences for nature and people [73,74,75,76]; it is important to keep monitoring forest EVI changes in the AGPE to better understand the confounding effects of human activities on forests.

The relatively low importance of forest composition in our analysis was unexpected as this variable is key for understanding vegetation and carbon dynamics [34,77,78]. The effect of this variable was likely lessened in the modeling given that the values represent means within a 250 m pixel rather than values for single trees. Similarly, the effect of forest protection status was not that important for explaining break magnitude variation compared to climate variables and only marginally important for positive breaks compared to negative ones. This suggests that most breaks are related to broad-acting natural disturbances such as climate stress that are likely to affect most forests within a region regardless of protection status.

Monitoring the state of forests with earth observation (EO) data is continually changing and expanding. There are data from different remote sensing sensors as well as a multitude of remote-sensing derived indices, each with different spatiotemporal resolutions, extents, and uses [24,26,27,79]. Similarly, there are several change detection algorithms in the literature each with its own pros and cons. For example, Militino et al. [62], Pasquarella et al. [80], and Yan et al. [81] compared several change detection algorithms including BFAST and LandTrendr. Here, we contribute to the science of EO forest monitoring by quantifying forest vegetation change, which can shed light on forest health and the capacity of forests to provide ecosystem services. We see the spatiotemporal analysis of abrupt and gradual change as an important tool in a growing forest monitoring toolbox that could include multi-algorithmic approaches [82] to improve change detection and attribution [19].

4.3. Limitations

A noteworthy limitation was the pixel resolution used (MODIS 250 m). This likely prevented the detection of disturbances that operate at finer spatial scales and affected matches between negative breaks and disturbance data whose original resolution was 30 m. Choosing a different h parameter would have produced more or fewer breaks, which in turn could have affected our matches with disturbance data and our interpretation of the results. Moreover, given our reliance on remote sensing data, our study focused on overstory forest vegetation dynamics with no attention to understory dynamics. Furthermore, EVI as well as other VIs, experience saturation when measured over dense canopies [67], and this may have prevented the detection of changes over time, especially in protected forest pixels. Lastly, the breaks and trends found here were not validated with forest plot data.

We have not linked our results directly to the provisioning of ecosystem services like carbon sequestration. This is because we did not assess the relationship between EVI breaks (or trends) with forest carbon data. However, it can be implied that greening forests are associated with higher biomass and forest ecosystem services compared to browning forests. Future work could relate EVI breaks and trends (positive and negative) to gross primary productivity (GPP) or to carbon loss/gain.

Our results should be interpreted with caution as they are derived from a single metric, EVI. A more complicated (but perhaps more robust) approach could have been to use other VIs in conjunction with EVI. Some studies have shown that NDVI and EVI can complement each other in certain landscapes [83]. Perhaps doing so would have improved the spatiotemporal matching with disturbance data. More importantly, here, we focused on a proxy of the ecological supply of forest ecosystem services. Monitoring other dimensions of forest ecosystem services such as anthropogenic contributions and instrumental value [21,22] could also reveal important changes in these services.

5. Conclusions

Forest EVI changes may result from processes that occurred in the recent past or spatially nearby (i.e., temporal and spatial lagged effects) and/or due to historical spatial legacies [84]. In this sense, EVI breaks and trends are not only due to multiple drivers but their relative contribution varies through time. While EVI mostly reflects forest canopy conditions, we found that most of the changes detected were negative breaks occurring mostly in non-protected areas and that not all negative breaks detected could be linked to disturbance data (i.e., forest fires, harvest, and insect damage). The negative and positive breaks detected were related to forest age and climate conditions (summer droughts and warm winters) and climate variables were found to have multi-year lag effects on EVI values. Protected forest pixels showed proportionally fewer negative and more positive breaks compared to non-protected forest pixels. Furthermore, on average, the magnitude of abrupt greening was higher in protected forest pixels compared to non-protected forest pixels suggesting that protected forests are healthier. The most common trend class in the study area involved greening with the highest proportion occurring in non-protected forests. Relatively large clusters of long-term browning trends were located in the northeast region of the study area. Further assessment of the spatiotemporal patterns of EVI breaks and trends is needed to shed light on the local drivers of forest health and their potential effects on forest ecosystem services and dynamics.

Author Contributions

Conceptualization: P.S.R., A.M.S. and M.-J.F.; methodology: P.S.R., A.M.S. and M.-J.F.; formal analysis, P.S.R.; data curation, P.S.R.; writing: P.S.R., M.-J.F., A.M.S. and A.G.; funding acquisition: M.-J.F. and A.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) for ResNet; (funding reference number NSERC NETGP 523374-18). Cette recherche a été financée par le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG); (numéro de référence NSERC NETGP 523374-18). Marie-Josée Fortin acknowledges funding from the NSERC of Canada Discovery Grant (no. 5134) and the NSERC Canada Research Chair (CRC). Andrew Gonzalez acknowledges support from the Liber Ero Chair in Biodiversity Conservation and the NSERC Discovery program.

Data Availability Statement

The original data presented in this study are openly available in the Land Processes Distributed Active Archive Center (LP DAAC) (doi:10.5067/MODIS/MOD13Q1.061).

Acknowledgments

We thank Jack Goldman for his valuable advice and suggestions. We also thank Flavio Affinito and Carina Rauen Firkowski for their support and feedback. Lastly, we thank the Digital Research Alliance of Canada for providing high-performance computing resources for our research.

Conflicts of Interest

Author Amanda M. Schwantes was employed by the company Apex Resource Management Solutions. The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

AGPE	Algonquin Park and Greater Park Ecosystem (Algonquin Provincial Park and
	surrounding area)
BFAST	Breaks For Additive Seasonal and Trend
BRT	boosted regression tree
CanLaD	Canada Landsat Disturbance data source
EVI	enhanced vegetation index
GFC	Global Forest Change
MODIS	Moderate Resolution Imaging Spectroradiometer
VI	vegetation index
XGBoost	extreme gradient boosting

Appendix A

Appendix A.1. Data Sources

Table A1. Main data sources used in this paper. Several data sources came from Canada’s National Terrestrial Ecosystem Monitoring System (NTEMS) https://opendata.nfis.org/mapserver/nfis-change_eng.html (accessed on 30 May 2023).

Data	Time Period	Spatial Resolution	Temporal Resolution	Source, Author
MODIS EVI	2003–2022	250 m	16-days	NASA Land Processes Distributed Active Archive Center (LP DAAC)
Forest fire	2003–2020	30 m	Yearly	Hermosilla et al. [46]
Forest harvest	2003–2020	30 m	Yearly	Hermosilla et al. [46]
Insect damage polygons	2003–2020	polygons	Yearly	Ontario Ministry of Natural Resources
CanLaD (fire + harvest)	2003–2015	30 m	Yearly	Guindon et al. [47]
Global Forest Change (GFC)	2003–2021	30 m	Yearly	Hansen et al. [48]
Climate variables from ClimateNA	2000–2020	10 km (downscaled to 250 m)	Yearly	Wang et al. [50]
Forest composition (percentage of conifers derived from these sources)	2001, 2011, 2019	250 m, 30 m	Year	Beaudoin et al. [43], Hermosilla et al. [42]
Forest age	2001, 2011, 2019	250 m, 30 m	Year	Beaudoin et al. [43], Maltman et al. [44]
Forest land cover	2003	30 m	Year	Hermosilla et al. [45]
Canada’s protected areas (polygons)	2020	polygons	-	Environment and Climate Change Canada
Elevation	-	30 m	-	NASA Shuttle Radar Topography Mission Global

The spatial resolution used for the response data (EVI) was 250 m. Raster data were either aggregated or downscaled to match the spatial resolution of EVI data. Polygon vector data were transformed to raster format to ease analysis. Forest age and percentage of conifers annual time series were estimated using pixel-based annual rates of change derived from the cross-sectional data available (i.e., 2001, 2011, and 2019). All data processing was conducted with R (version 4.1.1, [39]).

We extracted 30 climate-related variables from ClimateNA (version 7.31, [50]) for the period 1998–2021 using elevation raster centroids with a 250 m resolution. Climate variables were first extracted as text files from ClimateNA and then converted to raster time series using R (packages sqldf (version 0.4-11) and terra (version 1.6-17)). We applied the variance inflation factor (<5) to keep a smaller set of less correlated variables: winter degree days above 5 °C (winter DD5) and summer climate moisture index (summer CMI). We also used Canada’s protected areas boundaries (Environment and Climate Change Canada [2] and elevation data (NASA Shuttle Radar Topography Mission Global, ≈30 m resolution, product: SRTMGL1 v003).

Appendix A.2. Break and Trend Detection with BFAST Algorithms

EVI raster data were cleaned prior to running BFAST’s algorithms following the cleaning steps proposed by Samanta et al. [41]. Time series for the period between 2003 and 2022 were then derived for each pixel and used as inputs to bfast and bfast01 functions. The magnitude, type (or direction), and time of breaks (including 95% confidence intervals, CIs) were derived from bfast. The magnitude of breaks was used as the response variable in boosted regression tree models as described below. The bfast01 function provided inputs to the bfastclassify function, which classified each time series into one of eight trend types [57]. The outputs of bfast and bfastclassify were then converted back to raster files. A version of the R code used here can be found on GitHub (https://github.com/georod/forc_trends) (accessed on 24 June 2024).

Apart from the sensitivity analyses shown in Table A2 and Table A3, where the values of bfast’s h parameter varied, all other analyses presented here used h = 0.05. This value of h is equivalent to 23 observations or 1 year of data given the temporal resolution (16 days) and length (20 years, 2003–2022) of the EVI time series.

Table A2. Sensitivity analysis of bfast’s h parameter. This parameter determines the temporal window size (or minimal segment size) between potential breaks. The temporal equivalent of h is shown in brackets along the top row. For example, h = 0.05 is equivalent to one year of observations given the 16-day temporal resolution and length of the time series used here (20 years, 2003–2022). The number of breaks, for both positive and negative breaks, changed with different values of h. The highest number of breaks was produced when h was 0.05. However, the highest value of breaks per pixel (2.50) and the highest proportion of positive breaks (44%) were produced when h was 0.025 (6 months). A given value of h will be good for detecting some types of disturbance but not all (Fang et al. [60]). For the main analyses, h = 0.05 was chosen as this allowed joining breaks with yearly disturbance data and the production of annual summaries.

	h = 0.025	h = 0.05	h = 0.1	h = 0.15
	(6 Months)	(1 Year)	(2 Years)	(3 Years)
Total breaks	1477	15,904	11,799	10,701
Unique pixels	590	12,348	10,146	9662
Breaks per pixel	2.50	1.29	1.16	1.11
Negative breaks (%)	822 (56)	11,969 (75)	9372 (79)	8466 (79)
Unique pixels	440	9720	8276	7792
Breaks per pixel	1.87	1.23	1.13	1.09
Positive breaks (%)	655 (44)	3935 (25)	2427 (21)	2235 (21)
Unique pixels	308	3033	2115	2057
Breaks per pixel	2.13	1.30	1.15	1.09

Table A3. Matching success (overlap) between forest EVI breaks and various disturbance data sources. Matches across several values of h (a key BFAST parameter) are shown. The highest percentage of negative and positive matches occurred when h was 0.025 (6 months). The percentage matches were comparable for positive breaks when h was 0.05, 0.1, and 0.15 and higher for negative breaks for the latter two. Matches with positive breaks can be seen as an index of the false positive rate as disturbance data should not match with positive breaks (abrupt greenness). Matches were conducted based on pixel and time of break 95% confidence intervals overlap. The data sources listed are those shown in Appendix A.1 above. An h = 0.05 was used in the main analysis as this allowed joining breaks with yearly disturbance data and the production of annual summaries.

Break Type	Data Source	h = 0.05		h = 0.025		h = 0.1		h = 0.15
		(1 Year)		(6 Months)		(2 Years)		(3 Years)
		Breaks Matched	Percent Matched	Breaks Matched	Percent Matched	Breaks Matched	Percent Matched	Breaks Matched	Percent Matched
Negative
	Fire	148	1.2	44	5.4	145	1.5	143	1.7
	Harvest	1864	15.6	246	29.9	1773	18.9	1655	19.5
	Insects	54	0.5	1	0.1	55	0.6	52	0.6
	CanLaD	3358	28.1	336	40.9	3100	33.1	2909	34.4
	GFC	7228	60.4	541	65.8	6419	68.5	5981	70.6
Positive
	Fire	0	0.0	2	0.3	0	0.0	1	0.0
	Harvest	189	4.8	122	18.6	91	3.7	82	3.7
	Insects	49	1.2	0	0.0	31	1.3	32	1.4
	CanLaD	694	17.6	256	39.1	449	18.5	406	18.2
	GFC	1272	32.3	324	49.5	768	31.6	693	31.0

Figure A1. Schematic of break and trend types employed in this study. A time series (ts) can be characterized by the presence or absence of breaks and trends. Breaks represent abrupt changes in a ts whereas trends represent gradual changes. Here, we refer to three major groups of trends: monotonic (a,b), interruption (c), and reversal trends (d). Similarly, we refer to two break types: negative (red downward arrow) and positive (green upward arrow) breaks. Trends can be monotonic, either increasing (green slopes) or decreasing (orange slopes). The former are referred to as greening trends and the latter as browning trends. Monotonic trends can show no breaks (a) or show one or more breaks (negative or positive) but in concordance with the slope of the trend segments (b). Conversely, interruption (c) and reversal (d) trends are characterized by having a break type in discordance with the slope of the trend segments. Interruption trends can have two positive trend segments divided by a negative break and vice versa (two negative trend segments divided by a positive break). Reversal trends have opposite trend segments divided by a negative or positive break. Lastly, some ts may not change or show changes that are too small to be detected with the methods employed (horizontal gray dashed line in (a)). Here, we use the trend classification proposed by de Jong et al. [57]—MIG: monotonic increasing, greening trend (without breaks) (bottom line in (a)); MDB: monotonic decreasing, browning trend (without breaks) (top line in (a)); MIGpb: monotonic increasing, greening trend with a positive break (top set of lines in (b)); MDBnb: monotonic decreasing, browning trend with a negative break (bottom set of lines in (b)); IInb: interruption, the increasing trend with a negative break (top set of lines in (c)); IDpb: interruption, decreasing trend with a positive break (bottom set of lines in (c)); RGB: reversal, greening to browning trend (top two sets of lines in (d)); RBG: reversal, browning to greening trend (bottom two sets of lines in (d)).

Figure A2. Density plots of forest EVI magnitude of breaks in the AGPE from 2003 to 2022. The number of breaks and their magnitudes broken down by year are shown. Extreme magnitude values have been omitted to aid visualization. Vertical lines show medians—solid: yearly; dashed: entire time series. The total number of breaks was 15,904 (11,969 negative and 3935 positive). The time of break was rounded up prior to plotting which caused 2003 breaks (8 negative and 18 positive) to be part of 2004. No breaks were detected in 2022.

Appendix A.3. Boosted Regression Trees (XGBoost Models)

The data frame used for running boosted regression trees (BRT) consisted of EVI break magnitudes as the response variable and forest and climate-related variables as the predictors. We used variation inflation factor (VIF) to simplify a first set of yearly climate predictors (n = 30) and their corresponding temporal lagged versions (1-, 2-, and 3-year lags). Only climate variables with VIF < 5 were kept, removing highly collinear variables. The final set of predictors (n = 9) used in modeling included forest age, percentage of conifer composition, forest protection status, winter degree days above 5 °C (1- and 3-year lags), and summer climate moisture index (1-, 2-, and 3-year lags). Every row in the data frame was a break matched with predictors by pixel ID and (floored) year of break. The joining of response and predictor records was conducted with SQL using DuckDB (version 0.10.3).

For most breaks, the 95% confidence intervals (95% CI) for the timing of breaks ranged from 1 to >80 months. We ran the BRT models on subsets of the data reducing the noise and improving model fit. Separate BRT models were run for negative and positive EVI breaks. For negative breaks, a model that only included records with 95% CI range < 6 months and break magnitudes >|−500| (n = 116), resulted in lower RMSE and higher R² compared to a model that used all negative break records (Table A4). For positive breaks, a model that only included records with break magnitudes > 500 and that did overlap with disturbance data (GFC data) showed the highest R² (n = 3263). Filtering out positive records that had wide 95% CI ranges did not improve the positive break model.

BRT modeling was conducted using the Extreme Gradient Boosting (XGBoost) library (https://xgboost.readthedocs.io/en/stable/index.html (accessed on 29 July 2023)) and the scikit-learn API interface (https://xgboost.readthedocs.io/en/stable/python/python_api.html#module-xgboost.sklearn (accessed on 29 July 2023)). Key XGBoost parameters and their corresponding values are shown below (Table A5). All other model parameters used default values (https://xgboost.readthedocs.io/en/stable/parameter.html (accessed on 4 June 2024)). Values for learning rate, max depth, and gamma were fine-tuned using a randomized search approach (function RandomizedSearchCV in scikit-learn) to reduce over-fitting. Model evaluation was conducted using a train and test set approach, with the test having 33% of the records. XGBoost (version 1.7.3) was run with Python (version 3.10) and Jupyter Notebooks. The Python code used here can be found on GitHub (https://github.com/georod/forc_trends/ (accessed on 24 June 2024)).

Table A4. Boosted regression tree (BRT) model performance. Two separate models were run, one for negative breaks and another for positive breaks. "The number of records differed but the same variables were used in both models. For each type of break, two models are shown: all records and a subset of records. The former includes all the breaks detected while the latter includes a subset of all breaks. For negative breaks, the subset only included records with narrow time of break 95% CI ranges (<1.5 months) and magnitude of change >|−500| (n = 116). For positive breaks, the subset only included records with a magnitude of change >500 that did not overlap with GFC disturbance data (n = 3263). Models were evaluated using train (67% of records) and test (33% of records) sets. BRT was run using the XGBoost library the scikit-learn API interface.

Model	MSE	RMSE	R²	R²-adj	n Records	n Variables
Negative break magnitudes
All records	81,922.19	286.22	0.334	0.332	11,969	9
Subset of records	85,546.93	292.48	0.724	0.639	116	9
Positive break magnitudes
All records	100,565.73	317.12	0.337	0.333	3935	9
Subset of records	96,362.44	310.42	0.358	0.352	3263	9

Table A5. Key XGBoost parameters and their corresponding values as used in the modeling. The parameter values were the same for negative and positive break models. XGBoost was the library used to run boosted regression trees together with the scikit-learn API interface.

XGBoost Parameter	XGBoost Parameter Value
booster	gbtree
early_stopping_rounds	50
gamma	0.2
learning_rate	0.01
max_depth	8
n_estimators	1000
random_state	42
reg_lambda	10
reg_alpha	1

Figure A3. Ranking of all predictors (features) used in the XGBoost models. Panel (a) shows predictors of negative break magnitudes and panel (b) shows those of positive break magnitudes. These models were run with subsets of the detected breaks (negative breaks, n = 116 records; positive, n = 3263 records). Variable abbreviations—for_age: forest age; dd5_wt: winter degree days above 5 °C; cmi_sm: summer climate moisture index; for_con: percentage of conifers; for_pro_0: non-protected forest; lag#: 1-, 2- or 3-year lags. The protected forest variable is not present in (a) given its lack of importance in explaining negative breaks.

Figure A4. Partial dependence plots of magnitude of negative break predictors (features). The relationships between predictors and response (magnitude of EVI breaks) variables were mostly non-linear. Plots were created from the output of XGBoost. The model was run with a subset of all detected negative breaks (n = 116). The values on the y-axes are absolute values of negative magnitudes. Variable abbreviations—for_age: forest age; dd5_wt: winter degree days above 5 °C; cmi_sm: summer climate moisture index; for_con: percentage of conifers; for_pro_0: non-protected forest equals 1, protected equals 0; lag#: 1-, 2- or 3-year lags.

Figure A5. Partial dependence plots of magnitude of positive break predictors (features). The relationships between predictors and response (magnitude of EVI breaks) variables were mostly non-linear. Plots were created from the output of XGBoost. The model was run with a subset of all the detected positive breaks (n = 3263). Variable abbreviations—for_age: forest age; dd5_wt: winter degree days above 5 °C; cmi_sm: summer climate moisture index; for_con: percentage of conifers; for_pro_0: non-protected forest equals 1, protected equals 0; lag#: 1-, 2- or 3-year lags.

Appendix A.4. Geographic Distribution of Breaks and Trends

Figure A6. Geographic distribution of trends in Algonquin Park and the surrounding area, which we refer to as the Algonquin Greater Park Ecosystem (AGPE). All maps show trends that were derived from the output of bfastclassify. Compared to greening trends (MIG) which occurred throughout the AGPE (a), browning trends (MDB) mostly occurred in the NE quadrant (b). Most increasing trends with negative breaks (interruptions, IInb) occurred in the NW quadrant (c) while most of the relatively few decreasing trends with positive breaks (interruptions, IDpb) occurred in the NE quadrant (c). Notably, browning to greening reverse trends (RBG) co-occurred with browning trends in the NE quadrant (d). The dashed black line shows the study area footprint (1.5 M ha). The solid black line shows Algonquin Provincial Park’s footprint. The perpendicular dashed lines divide the study area into four quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the description of the results. Map projection across all figures is NAD83 Statistics Canada Lambert, EPSG: 3347.

Appendix A.5. Comparison of Disturbance Data

Table A6. Comparison of disturbance data sources used in this study. The percentages represent the amount of overlap (space and time) between data sources. When a data source is compared to itself it has a 100% overlap. As seen, there are discrepancies between different data sources. For instance, although CanLaD [47] captures fire and harvest disturbance pixels, only 17.2% of fire and 48.2% of harvest pixels from Hermosilla et al. [46] overlapped with it. GFC [48], which is a global non-labeled disturbance data source, overlapped with a higher amount of fire (40.5%), harvest (63.4%), and CanLaD (56.9%) pixels. Some level of discrepancy is expected, as the assumptions and methods used to derive each disturbance data source are different. The data sources are those shown in Table A1.

	Fire	Harvest	Insects	CanLaD	GFC
Fire	100.0%	0.2%	0.0%	0.2%	0.2%
Harvest	4.5%	100%	0.3%	14.8%	6.0%
Insects	0.0%	0.2%	100%	0.2%	0.4%
CanLaD	17.2%	48.2%	1.1%	100.0%	17.6%
GFC	40.5%	63.4%	7.3%	56.9%	100.0%

References

Hanes, C.C.; Wang, X.; Jain, P.; Parisien, M.A.; Little, J.M.; Flannigan, M.D. Fire-regime changes in Canada over the last half century. Can. J. For. Res. 2019, 49, 256–269. [Google Scholar] [CrossRef]
Canadian Forest Service. The State of Canada’s Forests: Annual Report 2021; Technical Report; Natural Resources Canada: Ottawa, ON, Canada, 2022. [Google Scholar]
Amiro, B.D.; Todd, J.B.; Wotton, B.M.; Logan, K.A.; Flannigan, M.D.; Stocks, B.J.; Mason, J.A.; Martell, D.L.; Hirsch, K.G. Direct carbon emissions from Canadian forest fires, 1959–1999. Can. J. For. Res. 2001, 31, 512–525. [Google Scholar] [CrossRef]
Kurz, W.A.; Dymond, C.C.; Stinson, G.; Rampley, G.J.; Neilson, E.T.; Carroll, A.L.; Ebata, T.; Safranyik, L. Mountain pine beetle and forest carbon feedback to climate change. Nature 2008, 452, 987–990. [Google Scholar] [CrossRef]
Hartmann, H.; Moura, C.F.; Anderegg, W.R.L.; Ruehr, N.K.; Salmon, Y.; Allen, C.D.; Arndt, S.K.; Breshears, D.D.; Davi, H.; Galbraith, D.; et al. Research frontiers for improving our understanding of drought-induced tree and forest mortality. New Phytol. 2018, 218, 15–28. [Google Scholar] [CrossRef]
Au, J.; Bloom, A.A.; Parazoo, N.C.; Deans, R.M.; Wong, C.Y.S.; Houlton, B.Z.; Magney, T.S. Forest productivity recovery or collapse? Model-data integration insights on drought-induced tipping points. Glob. Chang. Biol. 2023, 29, 5652–5665. [Google Scholar] [CrossRef] [PubMed]
Liu, Q.; Peng, C.; Schneider, R.; Cyr, D.; McDowell, N.G.; Kneeshaw, D. Drought-induced increase in tree mortality and corresponding decrease in the carbon sink capacity of Canada’s boreal forests from 1970 to 2020. Glob. Chang. Biol. 2023, 29, 2274–2285. [Google Scholar] [CrossRef]
Jactel, H.; Petit, J.; Desprez-Loustau, M.L.; Delzon, S.; Piou, D.; Battisti, A.; Koricheva, J. Drought effects on damage by forest insects and pathogens: A meta-analysis. Glob. Chang. Biol. 2012, 18, 267–276. [Google Scholar] [CrossRef]
Sommerfeld, A.; Senf, C.; Buma, B.; D’Amato, A.W.; Després, T.; Díaz-Hormazábal, I.; Fraver, S.; Frelich, L.E.; Gutiérrez, Á.G.; Hart, S.J.; et al. Patterns and drivers of recent disturbances across the temperate forest biome. Nat. Commun. 2018, 9, 4355. [Google Scholar] [CrossRef] [PubMed]
Seidl, R.; Schelhaas, M.J.; Rammer, W.; Verkerk, P.J. Increasing forest disturbances in Europe and their impact on carbon storage. Nat. Clim. Chang. 2014, 4, 806–810. [Google Scholar] [CrossRef]
Seidl, R.; Spies, T.A.; Peterson, D.L.; Stephens, S.L.; Hicke, J.A. Searching for resilience: Addressing the impacts of changing disturbance regimes on forest ecosystem services. J. Appl. Ecol. 2016, 53, 120–129. [Google Scholar] [CrossRef]
Bunn, A.G.; Goetz, S.J.; Kimball, J.S.; Zhang, K. Northern high-latitude ecosystems respond to climate change. EOS Trans. Am. Geophys. Union 2007, 88, 333–335. [Google Scholar] [CrossRef]
Périé, C.; de Blois, S. Dominant forest tree species are potentially vulnerable to climate change over large portions of their range even at high latitudes. PeerJ 2016, 4, e2218. [Google Scholar] [CrossRef] [PubMed]
Brice, M.H.; Cazelles, K.; Legendre, P.; Fortin, M.J. Disturbances amplify tree community responses to climate change in the temperate–boreal ecotone. Glob. Ecol. Biogeogr. 2019, 28, 1668–1681. [Google Scholar] [CrossRef]
Gonzalez, A.; Londoño, M.C. Monitor biodiversity for action. Science 2022, 378, 1147. [Google Scholar] [CrossRef] [PubMed]
Buxton, R.T.; Bennett, J.R.; Reid, A.J.; Shulman, C.; Cooke, S.J.; Francis, C.M.; Nyboer, E.A.; Pritchard, G.; Binley, A.D.; Avery-Gomm, S.; et al. Key information needs to move from knowledge to action for biodiversity conservation in Canada. Biol. Conserv. 2021, 256, 108983. [Google Scholar] [CrossRef]
Lindenmayer, D.B.; Likens, G.E. Adaptive monitoring: A new paradigm for long-term research and monitoring. Trends Ecol. Evol. 2009, 24, 482–486. [Google Scholar] [CrossRef]
McDowell, N.G.; Coops, N.C.; Beck, P.S.A.; Chambers, J.Q.; Gangodagamage, C.; Hicke, J.A.; Huang, C.Y.; Kennedy, R.; Krofcheck, D.J.; Litvak, M.; et al. Global satellite monitoring of climate-induced vegetation disturbances. Trends Plant Sci. 2015, 20, 114–123. [Google Scholar] [CrossRef]
Gonzalez, A.; Chase, J.M.; O’Connor, M.I. A framework for the detection and attribution of biodiversity change. Philos. Trans. R. Soc. B Biol. Sci. 2023, 378, 20220182. [Google Scholar] [CrossRef]
Kissling, W.D.; Hardisty, A.; García, E.A.; Santamaria, M.; De Leo, F.; Pesole, G.; Freyhof, J.; Manset, D.; Wissel, S.; Konijn, J.; et al. Towards global interoperability for supporting biodiversity research on essential biodiversity variables (EBVs). Biodiversity 2015, 16, 99–107. [Google Scholar] [CrossRef]
Balvanera, P.; Brauman, K.A.; Cord, A.F.; Drakou, E.G.; Geijzendorffer, I.R.; Karp, D.S.; Martín-López, B.; Mwampamba, T.H.; Schröter, M. Essential ecosystem service variables for monitoring progress towards sustainability. Curr. Opin. Environ. Sustain. 2022, 54, 101152. [Google Scholar] [CrossRef]
Schwantes, A.M.; Firkowski, C.R.; Affinito, F.; Rodriguez, P.S.; Fortin, M.J.; Gonzalez, A. Monitoring ecosystem services with Essential Ecosystem Service Variables. Front. Ecol. Environ. 2024. [Google Scholar] [CrossRef]
Kerr, J.T.; Ostrovsky, M. From space to species: Ecological applications for remote sensing. Trends Ecol. Evol. 2003, 18, 299–305. [Google Scholar] [CrossRef]
Xue, J.; Su, B. Significant Remote Sensing Vegetation Indices: A Review of Developments and Applications. J. Sens. 2017, 2017, 1–17. [Google Scholar] [CrossRef]
Lechner, A.M.; Foody, G.M.; Boyd, D.S. Applications in Remote Sensing to Forest Ecology and Management. ONE Earth 2020, 2, 405–412. [Google Scholar] [CrossRef]
Huete, A.R. Vegetation Indices, Remote Sensing and Forest Monitoring. Geogr. Compass 2012, 6, 513–532. [Google Scholar] [CrossRef]
Petropoulos, G.P.; Kalaitzidis, C. Multispectral vegetation indices in remote sensing: An overview. In Ecological Modeling; Zhang, W.J., Ed.; Environmental Science, Engineering and Technology, Nova Science Publishers: Hauppauge, NY, USA, 2012; Chapter 2; pp. 15–39. [Google Scholar]
Beaudoin, A.; Bernier, P.Y.; Guindon, L.; Villemaire, P.; Guo, X.J.; Stinson, G.; Bergeron, T.; Magnussen, S.; Hall, R.J. Mapping attributes of Canada’s forests at moderate resolution through kNN and MODIS imagery. Can. J. For. Res. 2014, 44, 521–532. [Google Scholar] [CrossRef]
Cartwright, J.M.; Littlefield, C.E.; Michalak, J.L.; Lawler, J.J.; Dobrowski, S.Z. Topographic, soil, and climate drivers of drought sensitivity in forests and shrublands of the Pacific Northwest, USA. Sci. Rep. 2020, 10, 18486. [Google Scholar] [CrossRef] [PubMed]
Wu, Z.; Liu, Y.; Han, Y.; Zhou, J.; Liu, J.; Wu, J. Mapping farmland soil organic carbon density in plains with combined cropping system extracted from NDVI time-series data. Sci. Total Environ. 2021, 754, 142120. [Google Scholar] [CrossRef] [PubMed]
Sothe, C.; Gonsamo, A.; Arabian, J.; Kurz, W.A.; Finkelstein, S.A.; Snider, J. Large Soil Carbon Storage in Terrestrial Ecosystems of Canada. Glob. Biogeochem. Cycles 2022, 36, e2021GB007213. [Google Scholar] [CrossRef]
Jentsch, A.; White, P. A theory of pulse dynamics and disturbance in ecology. Ecology 2019, 100, e02734. [Google Scholar] [CrossRef]
Coops, N.C.; Shang, C.; Wulder, M.A.; White, J.C.; Hermosilla, T. Change in forest condition: Characterizing non-stand replacing disturbances using time series satellite imagery. For. Ecol. Manag. 2020, 474, 118370. [Google Scholar] [CrossRef]
Silva Pedro, M.; Rammer, W.; Seidl, R. Tree species diversity mitigates disturbance impacts on the forest carbon cycle. Oecologia 2015, 177, 619–630. [Google Scholar] [CrossRef] [PubMed]
Arain, M.A.; Xu, B.; Brodeur, J.J.; Khomik, M.; Peichl, M.; Beamesderfer, E.; Restrepo-Couple, N.; Thorne, R. Heat and drought impact on carbon exchange in an age-sequence of temperate pine forests. Ecol. Process. 2022, 11, 7. [Google Scholar] [CrossRef] [PubMed]
Ontario Parks. Algonquin Provincial Park Management Plan; Technical Report; Harpell Printing Ottawa Inc.: Ottawa, ON, Canada, 1998. [Google Scholar]
Strickland, D. Trees of Algonquin Park; Technical Report; The Friends of Algonquin Park: Whitney, ON, Canada, 1996. [Google Scholar]
Kremen, C.; Merenlender, A.M. Landscapes that work for biodiversity and people. Science 2018, 362, eaau6020. [Google Scholar] [CrossRef]
R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021. [Google Scholar]
Busetto, L.; Ranghetti, L. MODIStsp: An R package for automatic preprocessing of MODIS Land Products time series. Comput. Geosci. 2016, 97, 40–48. [Google Scholar] [CrossRef]
Samanta, A.; Ganguly, S.; Hashimoto, H.; Devadiga, S.; Vermote, E.; Knyazikhin, Y.; Nemani, R.R.; Myneni, R.B. Amazon forests did not green-up during the 2005 drought. Geophys. Res. Lett. 2010, 37. [Google Scholar] [CrossRef]
Hermosilla, T.; Wulder, M.A.; White, J.C.; Coops, N.C. Land cover classification in an era of big and open data: Optimizing localized implementation and training data selection to improve mapping outcomes. Remote Sens. Environ. 2022, 268, 112780. [Google Scholar] [CrossRef]
Beaudoin, A.; Bernier, P.Y.; Villemaire, P.; Guindon, L.; Guo, X.J. Tracking forest attributes across Canada between 2001 and 2011 using a k nearest neighbors mapping approach applied to MODIS imagery. Can. J. For. Res. 2018, 48, 85–93. [Google Scholar] [CrossRef]
Maltman, J.C.; Hermosilla, T.; Wulder, M.A.; Coops, N.C.; White, J.C. Estimating and mapping forest age across Canada’s forested ecosystems. Remote Sens. Environ. 2023, 290, 113529. [Google Scholar] [CrossRef]
Hermosilla, T.; Bastyr, A.; Coops, N.C.; White, J.C.; Wulder, M.A. Mapping the presence and distribution of tree species in Canada’s forested ecosystems. Remote Sens. Environ. 2022, 282, 113276. [Google Scholar] [CrossRef]
Hermosilla, T.; Wulder, M.A.; White, J.C.; Coops, N.C.; Hobart, G.W.; Campbell, L.B. Mass data processing of time series Landsat imagery: Pixels to data products for forest monitoring. Int. J. Digit. Earth 2016, 9, 1035–1054. [Google Scholar] [CrossRef]
Guindon, L.; Bernier, P.; Gauthier, S.; Stinson, G.; Villemaire, P.; Beaudoin, A. Missing forest cover gains in boreal forests explained. Ecosphere 2018, 9, e02094. [Google Scholar] [CrossRef]
Hansen, M.C.; Potapov, P.V.; Moore, R.; Hancher, M.; Turubanova, S.A.; Tyukavina, A.; Thau, D.; Stehman, S.V.; Goetz, S.J.; Loveland, T.R.; et al. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science 2013, 342, 850–853. [Google Scholar] [CrossRef]
Ontario Ministry of Natural Resources and Forestry. Forest Insect Damage Event. 2020. Available online: https://geohub.lio.gov.on.ca/documents/lio::forest-insect-damage-event/about (accessed on 3 March 2023).
Wang, T.; Hamann, A.; Spittlehouse, D.; Carroll, C. Locally Downscaled and Spatially Customizable Climate Data for Historical and Future Periods for North America. PLoS ONE 2016, 11, e0156720. [Google Scholar] [CrossRef] [PubMed]
Allen, C.D.; Macalady, A.K.; Chenchouni, H.; Bachelet, D.; McDowell, N.; Vennetier, M.; Kitzberger, T.; Rigling, A.; Breshears, D.D.; Hogg, E.T.; et al. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For. Ecol. Manag. 2010, 259, 660–684. [Google Scholar] [CrossRef]
Allen, C.D.; Breshears, D.D.; McDowell, N.G. On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere 2015, 6, 1–55. [Google Scholar] [CrossRef]
Environment and Climate Change Canada (ECCC). Canadian Protected and Conserved Areas Database (CPCAD)-Protected and Conserved Area Database. 2022. Available online: https://www.canada.ca/en/environment-climate-change/services/national-wildlife-areas/protected-conserved-areas-database.html (accessed on 25 May 2023).
Verbesselt, J.; Hyndman, R.; Newnham, G.; Culvenor, D. Detecting trend and seasonal changes in satellite image time series. Remote Sens. Environ. 2010, 114, 106–115. [Google Scholar] [CrossRef]
Verbesselt, J.; Zeileis, A.; Herold, M. Near real-time disturbance detection using satellite image time series. Remote Sens. Environ. 2012, 123, 98–108. [Google Scholar] [CrossRef]
de Jong, R.; de Bruin, S. Linear trends in seasonal vegetation time series and the modifiable temporal unit problem. Biogeosciences 2012, 9, 71–77. [Google Scholar] [CrossRef]
de Jong, R.; Verbesselt, J.; Zeileis, A.; Schaepman, M.E. Shifts in Global Vegetation Activity Trends. Remote Sens. 2013, 5, 1117–1133. [Google Scholar] [CrossRef]
Horion, S.; Prishchepov, A.V.; Verbesselt, J.; de Beurs, K.; Tagesson, T.; Fensholt, R. Revealing turning points in ecosystem functioning over the Northern Eurasian agricultural frontier. Glob. Chang. Biol. 2016, 22, 2801–2817. [Google Scholar] [CrossRef] [PubMed]
Lu, M.; Pebesma, E.; Sanchez, A.; Verbesselt, J. Spatio-temporal change detection from multidimensional arrays: Detecting deforestation from MODIS time series. ISPRS J. Photogramm. Remote Sens. 2016, 117, 227–236. [Google Scholar] [CrossRef]
Fang, X.; Zhu, Q.; Ren, L.; Chen, H.; Wang, K.; Peng, C. Large-scale detection of vegetation dynamics and their potential drivers using MODIS images and BFAST: A case study in Quebec, Canada. Remote Sens. Environ. 2018, 206, 391–402. [Google Scholar] [CrossRef]
Tong, X.; Brandt, M.; Yue, Y.; Horion, S.; Wang, K.; Keersmaecker, W.D.; Tian, F.; Schurgers, G.; Xiao, X.; Luo, Y.; et al. Increased vegetation growth and carbon stock in China karst via ecological engineering. Nat. Sustain. 2018, 1, 44–50. [Google Scholar] [CrossRef]
Militino, A.F.; Moradi, M.; Ugarte, M.D. On the Performances of Trend and Change-Point Detection Methods for Remote Sensing Data. Remote Sens. 2020, 12, 1008. [Google Scholar] [CrossRef]
de Jong, R.; Verbesselt, J.; Schaepman, M.E.; de Bruin, S. Trend changes in global greening and browning: Contribution of short-term trends to longer-term change. Glob. Chang. Biol. 2011, 18, 642–655. [Google Scholar] [CrossRef]
Elith, J.; Leathwick, J.R.; Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 2008, 77, 802–813. [Google Scholar] [CrossRef] [PubMed]
Chen, T.; Guestrin, C. XGBoost: A Scalable Tree Boosting System. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, 13–17 August 2016. [Google Scholar] [CrossRef]
Sharma, T.; Kurz, W.A.; Stinson, G.; Pellatt, M.G.; Li, Q. A 100-year conservation experiment: Impacts on forest carbon stocks and fluxes. For. Ecol. Manag. 2013, 310, 242–255. [Google Scholar] [CrossRef]
Gao, S.; Zhong, R.; Yan, K.; Ma, X.; Chen, X.; Pu, J.; Gao, S.; Qi, J.; Yin, G.; Myneni, R.B. Evaluating the saturation effect of vegetation indices in forests using 3D radiative transfer simulations and satellite observations. Remote Sens. Environ. 2023, 295, 113665. [Google Scholar] [CrossRef]
Bolton, D.K.; Coops, N.C.; Hermosilla, T.; Wulder, M.A.; White, J.C.; Ferster, C.J. Uncovering regional variability in disturbance trends between parks and greater park ecosystems across Canada (1985-2015). Sci. Rep. 2019, 9, 1323. [Google Scholar] [CrossRef]
Ghazoul, J.; Burivalova, Z.; Garcia-Ulloa, J.; King, L.A. Conceptualizing forest degradation. Trends Ecol. Evol. 2015, 30, 622–632. [Google Scholar] [CrossRef] [PubMed]
Stephens, B.B.; Gurney, K.R.; Tans, P.P.; Sweeney, C.; Peters, W.; Bruhwiler, L.; Ciais, P.; Ramonet, M.; Bousquet, P.; Nakazawa, T.; et al. Weak Northern and Strong Tropical Land Carbon Uptake from Vertical Profiles of Atmospheric CO₂. Science 2007, 316, 1732–1735. [Google Scholar] [CrossRef] [PubMed]
Michaelian, M.; Hogg, E.H.; Hall, R.J.; Arsenault, E. Massive mortality of aspen following severe drought along the southern edge of the Canadian boreal forest. Glob. Chang. Biol. 2011, 17, 2084–2094. [Google Scholar] [CrossRef]
Hartmann, H.; Schuldt, B.; Sanders, T.G.M.; Macinnis-Ng, C.; Boehmer, H.J.; Allen, C.D.; Bolte, A.; Crowther, T.W.; Hansen, M.C.; Medlyn, B.E.; et al. Monitoring global tree mortality patterns and trends. Report from the VW symposium ‘Crossing scales and disciplines to identify global trends of tree mortality as indicators of forest health’. New Phytol. 2018, 217, 984–987. [Google Scholar] [CrossRef]
Chiang, F.; Mazdiyasni, O.; AghaKouchak, A. Evidence of anthropogenic impacts on global drought frequency, duration, and intensity. Nat. Commun. 2021, 12, 2754. [Google Scholar] [CrossRef]
Hermans, K.; McLeman, R. Climate change, drought, land degradation and migration: Exploring the linkages. Curr. Opin. Environ. Sustain. 2021, 50, 236–244. [Google Scholar] [CrossRef]
Walker, D.W.; Van Loon, A.F. Droughts are coming on faster. Science 2023, 380, 130–132. [Google Scholar] [CrossRef]
Yuan, X.; Wang, Y.; Ji, P.; Wu, P.; Sheffield, J.; Otkin, J.A. A global transition to flash droughts under climate change. Science 2023, 380, 187–191. [Google Scholar] [CrossRef]
Mack, M.C.; Walker, X.J.; Johnstone, J.F.; Alexander, H.D.; Melvin, A.M.; Jean, M.; Miller, S.N. Carbon loss from boreal forest wildfires offset by increased dominance of deciduous trees. Science 2021, 372, 280–283. [Google Scholar] [CrossRef]
Augusto, L.; Boča, A. Tree functional traits, forest biomass, and tree species diversity interact with site properties to drive forest soil carbon. Nat. Commun. 2022, 13, 1097. [Google Scholar] [CrossRef]
Fassnacht, F.E.; White, J.C.; Wulder, M.A.; Næsset, E. Remote sensing in forestry: Current challenges, considerations and directions. For. Int. J. For. Res. 2023, 97, 11–37. [Google Scholar] [CrossRef]
Pasquarella, V.J.; Arévalo, P.; Bratley, K.H.; Bullock, E.L.; Gorelick, N.; Yang, Z.; Kennedy, R.E. Demystifying LandTrendr and CCDC temporal segmentation. Int. J. Appl. Earth Obs. Geoinf. 2022, 110, 102806. [Google Scholar] [CrossRef]
Yan, J.; He, H.; Wang, L.; Zhang, H.; Liang, D.; Zhang, J. Inter-Comparison of Four Models for Detecting Forest Fire Disturbance from MOD13A2 Time Series. Remote Sens. 2022, 14, 1446. [Google Scholar] [CrossRef]
Saxena, R.; Watson, L.T.; Wynne, R.H.; Brooks, E.B.; Thomas, V.A.; Zhiqiang, Y.; Kennedy, R.E. Towards a polyalgorithm for land use change detection. ISPRS J. Photogramm. Remote Sens. 2018, 144, 217–234. [Google Scholar] [CrossRef]
Huang, C.; Zhou, Z.; Wang, D.; Dian, Y. Monitoring forest dynamics with multi-scale and time series imagery. Environ. Monit. Assess. 2016, 188, 1–15. [Google Scholar] [CrossRef]
James, P.M.A.; Fortin, M.J.; Fall, A.; Kneeshaw, D.; Messier, C. The Effects of Spatial Legacies following Shifting Management Practices and Fire on Boreal Forest Age Structure. Ecosystems 2007, 10, 1261–1277. [Google Scholar] [CrossRef]

Figure 1. The study area is the Algonquin Provincial Park and the surrounding area, which we refer to as the Algonquin Greater Park Ecosystem (AGPE). Elevation (in meters above sea level, masl) is shown in (a). The forest’s mean age in 2019 is shown in (b). Protected areas within the AGPE are shown in (c). About 16% of forest pixels belong to protected areas. The geographic distribution of three disturbance agents in the period 2002–2020 is shown in (d). The gray color represents forested pixels and white non-forest pixels. The dashed black line shows the study area footprint (≈15,000 km² or 1.5 M ha). The solid black line shows Algonquin Provincial Park’s footprint. The two perpendicular dashed lines in (a) divide the study area into quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the description of the results. Map projection is NAD83 Statistics Canada Lambert, EPSG: 3347.

Figure 2. Main methodological steps used in this study. After downloading MODIS EVI data, pixels were filtered to only keep pixels with good quality data, within forested areas. EVI time series were then created and processed with three BFAST algorithms: bfast, bfast01, and bfastclassify. The bfast algorithm decomposes a time series into a seasonal component, trend, and noise. Using piece-wise linear regression on the trend component, it detects one or more breaks (if any). Here, we use three main outputs provided by bfast: type of break, magnitude of break, and time of break with 95% confidence intervals (CIs) The bfast01 algorithm runs a seasonally adjusted regression model on the ts and only detects the major break (if any). The bfastclassify algorithm then uses bfast01’s output to classify trends into one of eight possible trend types (Figure A1). Only the magnitude of break values estimated by bfast was used in boosted regression trees (XGBoost models) to explore their relationship with predictor variables (dashed box at bottom, Equation (1) in main text). Satellite icon from flaticon.com (accessed on 29 April 2022).

Figure 3. Spatial distribution of EVI negative and positive breaks in the AGPE from 2003 to 2022. Most of the breaks were found in the eastern half of the AGPE and more particularly in the northeast quadrant. There were 11,871 pixels with negative breaks (red pixels) and 3893 pixels with positive breaks (cyan pixels). These breaks were estimated with the bfast algorithm. The dashed black line shows the study area footprint. The solid black line shows Algonquin Provincial Park’s footprint. The two perpendicular dashed lines divide the study area into four quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the interpretation of results. Map projection is NAD83 Statistics Canada Lambert, EPSG: 3347.

Figure 4. Spatial distribution of EVI trend types in the AGPE from 2003 to 2022. These trends were produced with the bfastclassify algorithm. The trend types are those proposed by de Jong et al. [57]. Abbreviations—MIG: monotonic increasing, greening trend (n = 33,683); MDB: monotonic decreasing, browning trend (n = 4981); IInb: interruption, increasing trend with a negative break (n = 11,637); RBG: reversal, browning to greening trend (n = 11,654). (Trends MIGpb, MDBnb, IDpb, and RGB are not shown given their low percentages). The dashed black line shows the study area footprint. The solid black line shows Algonquin Provincial Park’s footprint. The two perpendicular dashed lines divide the study area into four quadrants (northwest, NW; northeast, NE; southwest, SW; and southeast, SE) to ease the interpretation of results. Map projection across all figures is NAD83 Statistics Canada Lambert, EPSG: 3347.

Figure 5. Predictors of the magnitude of EVI breaks (negative breaks in red and positive breaks in cyan tone boxes) from 2003 to 2022. The ranking reflects feature importance using the gain metric as estimated by XGBoost models. The same five predictors are in the top five but with slightly different rankings (connecting lines with slopes) except for the summer climate moisture index with a 3-year lag, which ranks fourth in both (connecting line with no slope). Forest protection status is low-ranking for both types of breaks. The XGBoost models were run with subsets of the detected breaks (negative breaks, n = 116 records; positive breaks, n = 3263 records).

Table 1. List of key variables and terms used in this study. For the two climate-related predictors (winter DD5 and summer CMI), 1-, 2- and 3-year lagged versions were created and used in boosted regression trees (BRTs).

Description	Short Form	Usage
Enhanced vegetation index	EVI	Input to BFAST functions (bfast & bfast01).
EVI abrupt change or break	break	Output of bfast function. Used in descriptive statistics.
EVI magnitude of break	magnitude	Output of bfast function. Figure A2. Response variable in BRT models. BRT variable name: magnitude.
EVI type of break	break type	Output of bfast function. Negative and positive breaks are produced. Used in descriptive statistics. Table 2, Figure 3 and Figure A2.
EVI time of break	time of break	Output of bfast function. This variable is associated with 95% CIs. Used in descriptive statistics and for joining breaks with other data sources
EVI gradual change or trend	trend	Output of bfastclassify function. Used in descriptive statistics. Trend classification proposed by de Jong et al. [57].
EVI monotonic increasing, greening trend (without breaks)	MIG	Output of bfastclassify function. Used in descriptive statistics. Table 3, Figure 4 and Figure A4.
EVI monotonic decreasing, browning trend (without breaks)	MDB	Output of bfastclassify function. Used in descriptive statistics. Table 3, Figure 4 and Figure A4.
EVI monotonic increasing, greening trend with a positive break	MIGpb	Output of bfastclassify function. Used in descriptive statistics. Table 3.
EVI monotonic increasing, browning trend with a negative break	MDBnb	Output of bfastclassify function. Used in descriptive statistics. Table 3.
EVI interruption, increasing trend with a negative break	IInb	Output of bfastclassify function. Used in descriptive statistics. Table 3, Figure 4 and Figure A4.
EVI interruption, decreasing trend with a positive break	IDpb	Output of bfastclassify function. Used in descriptive statistics. Table 3, Figure A4.
EVI reversal, greening to browning trend	RGB	Output of bfastclassify function. Used in descriptive statistics. Table 3.
EVI reversal, browning to greening trend.	RBG	Output of bfastclassify function. Used in descriptive statistics. Table 3, Figure 4 and Figure A4.
Forest age	age	Forest-related predictor used in BRT models. BRT variable name: for_age
Percent coniferous composition	percentage of conifers	Forest-related predictor used in BRT models. BRT variable name: for_con.
Forest protection status	protection status	Forest-related predictor used in BRT models. BRT variable name: for_pro_0.
Winter degree days above 5 °C	winter DD5	Climate-related predictor used in BRT models. Temporal lagged versions were used. BRT variable names: dd5_wt_lag1, dd5_wt_lag3.
Summer climate moisture index	summer CMI	Climate-related predictor used in BRT models. Temporal lagged versions were used. BRT variable names: cmi_sm, cmi_sm_lag1, cmi_sm_lag2, cmi_sm_lag3.

Table 2. Summary of EVI breaks (abrupt change) from 2003 to 2022 broken down by break type and forest protection status. The values were derived from the output of the bfast algorithm. Percentages in brackets are derived from the overall number of breaks while percentages in square brackets are derived from the overall number of forest pixels. A pixel can have one or more breaks. Therefore, the values in the negative and positive columns (last two columns) do not add up to the values in the ‘All’ column (third last column).

Group	# of Forest Pixels	# of Breaks in 2003–2022			# of Pixels with Breaks in 2003–2022
		(% of Breaks)			[% of Forest Pixels]
		All	Negative	Positive	All	Negative	Positive
AGPE	230,806	15,904	11,969	3935	12,348	9720	3033
			(75.3%)	(24.7%)	[5.3%]	[4.2%]	[1.3%]
Non-protected forests	194,046	15,178	11,499	3679	11,805	9366	2840
			(75.8%)	(24.2%)	[6.1%]	[4.8%]	[1.5%]
Protected forests	36,760	726	470	256	543	354	193
			(64.7%)	(35.3%)	[1.5%]	[1.0%]	[0.5%]

Table 3. Summary of EVI trend (gradual change) pixels from 2003 to 2022 broken down by trend type and forest protection status. The values were derived from the output of the bfastclassify algorithm. Percent of forest pixels is the number of forest pixels in each trend type divided by the overall number of forest pixels. The percentage of all trends is the number of trend pixels in each trend type divided by the overall number of trend pixels. Trend types are those proposed by de Jong et al. [57]. Abbreviations—MIG: monotonic increasing, greening trend (without breaks); MDB: monotonic decreasing, browning trend (without breaks); MIGpb: monotonic increasing, greening trend with a positive break; MDBnb: monotonic decreasing, browning trend with a negative break; IInb: interruption, increasing trend with a negative break; IDpb: interruption, decreasing trend with a positive break; RGB: reversal, greening to browning trend; RBG: reversal, browning to greening trend.

Group	# of Forest Pixels	# of Trend Pixels	Trend Types
			MIG	MDB	MIGpb	MDBnb	IInb	IDpb	RGB	RBG
AGPE	230,806	64,188	33,683	4981	73	17	11,637	1674	469	11,654
% of forest pixels		27.8	14.6	2.2	0.0	0.0	5.0	0.7	0.2	5.0
% of all trends			52.5	7.8	0.1	0.0	18.1	2.6	0.7	18.2
Non-protected forests	194,046	58,457	30,370	4366	69	15	10,691	1445	404	11,097
% of forest pixels		28.9	15.0	2.2	0.0	0.0	5.3	0.7	0.2	5.5
% of all trends			52.0	7.5	0.1	0.0	18.3	2.5	0.7	19.0
Protected forests	36,760	5731	3313	615	4	2	946	229	65	557
% of forest pixels		14.6	8.4	1.6	0.0	0.0	2.4	0.6	0.2	1.4
% of all trends			57.8	10.7	0.1	0.0	16.5	4.0	1.1	9.7

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Rodriguez, P.S.; Schwantes, A.M.; Gonzalez, A.; Fortin, M.-J. Monitoring Changes in the Enhanced Vegetation Index to Inform the Management of Forests. Remote Sens. 2024, 16, 2919. https://doi.org/10.3390/rs16162919

AMA Style

Rodriguez PS, Schwantes AM, Gonzalez A, Fortin M-J. Monitoring Changes in the Enhanced Vegetation Index to Inform the Management of Forests. Remote Sensing. 2024; 16(16):2919. https://doi.org/10.3390/rs16162919

Chicago/Turabian Style

Rodriguez, Peter S., Amanda M. Schwantes, Andrew Gonzalez, and Marie-Josée Fortin. 2024. "Monitoring Changes in the Enhanced Vegetation Index to Inform the Management of Forests" Remote Sensing 16, no. 16: 2919. https://doi.org/10.3390/rs16162919

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu