Characterizing over Four Decades of Forest Disturbance in Minnesota, USA

: Spatial information about disturbance driven patterns of forest structure and ages across landscapes provide a valuable resource for all land management e ﬀ orts including cross-ownership collaborative forest treatments and restoration. While disturbance events in general are known to impact stand characteristics, the agent of change may also inﬂuence recovery and the supply of ecosystem services. Our study utilizes the full extent of the Landsat archive to identify the timing, extent, magnitude, and agent, of the most recent fast disturbance event for all forested lands within Minnesota, USA. To account for the di ﬀ erences in the Landsat sensors through time, speciﬁcally the coarser spatial, spectral, and radiometric resolutions of the early MSS sensors, we employed a two-step approach, ﬁrst harmonizing spectral indices across the Landsat sensors, then applying a segmentation algorithm to ﬁt temporal trends to the time series to identify abrupt forest disturbance events. We further incorporated spectral, topographic, and land protection information in our classiﬁcation of the agent of change for all disturbance patches. After allowing two years for the time series to stabilize, we were able to identify the most recent fast disturbance events across Minnesota from 1974–2018 with a change versus no-change validation accuracy of 97.2% ± 1.9%, and higher omission (14.9% ± 9.3%) than commission errors (1.6% ± 1.9%) for the identiﬁcation of change patches. Our classiﬁcation of the agent of change exhibited an overall accuracy of 96.5% ± 1.9% with classes including non-disturbed forest, land conversion, ﬁre, ﬂooding, harvest, wind / weather, and other rare natural events. Individual class errors varied, but all class user and producer accuracies were above 78%. The unmatched nature of the Landsat archive for providing comparable forest attribute and change information across more than four decades highlights the value of the totality of the Landsat program to the larger geospatial, ecological research, and forest management communities.


Introduction
Spatially explicit information about stand-level forest attributes across land ownerships are valuable for regional forest assessments, coordinated forest management, and policy development. Landscape-level patterns of forest structure and stand age are largely driven by disturbance regimes [1]. Characteristics of forest disturbance events, including the duration, severity, and causal agent, will affect the residual forest structure and potential recovery trajectories [2], both of which are vital to forest management planning and restoration actions. Additionally, the agent of disturbance (e.g., harvest, fire, insects) has varying impacts on forest ecosystem services such as soil productivity [3], water quality [4], support of biodiversity [5], and carbon sequestration [6]. Thus, understanding creation of annual cloud-free spectral indices from 1972 to present. Studies that have incorporated MSS imagery into a complete time series analyses have mostly focused on predicting current forest structure and biomass [20], stand structural stages and canopy cover through time [11,25], or disturbances within relatively small geographic areas [26]. The efficacy of using this complete time series for disturbance mapping and agent attributions across a large area, such as the forested region of an entire state, remains untested.
Quantification of disturbance and associated agent attribution across over four decades may provide valuable information for forest managers who are frequently tasked with managing their land bases to sustain multiple ecosystem services (e.g., timber production, water quality, wildlife habitat, carbon sequestration) across large spatial extents and temporal domains. The objectives of this study were to: (1) utilize the complete Landsat archive  to identify the timing, magnitude of change, and extent of the most recent fast disturbance events across a state-wide extent, utilizing Minnesota forested lands as a case study area; (2) classify all stand-altering change events by the major disturbance agents of the region and assess sources of classification errors; and (3) identify spatial and temporal patterns of disturbance agents and their magnitudes of change by ecological regions and land ownerships to inform all lands management efforts.

Study Area
We focused our case study on the forested lands of Minnesota, USA, which encompasses a variety of ownerships that are approximately evenly split between private and public entities ( Figure 1). Forests cover approximately 1/3 of the land area in Minnesota and are predominantly located in the northern-northeastern portion of the state (Figure 1). A variety of cover types occur in the region, with aspen-birch (Populus spp.-Betula spp.) and spruce-fir (Picea mariana-Abies balsamea) types dominant, and to a lesser extent pine (Pinus resinosa, Pinus strobus, Pinus banksiana), oak (Quercus spp.), and northern hardwoods [27]. Wetlands, including forested bogs, peatlands, and swamps, are found extensively throughout the forested region of the state. A variety of land uses occur in Minnesota forests including timber production, mining operations, and recreation. The four ecological provinces of Minnesota (Figure 1), which are defined by their climate, geologic history, soils, hydrology, and associated vegetation, include: Eastern Broadleaf Forest Province (4.9 million ha); Laurentian Mixed Forest Province (9.3 million ha); Prairie Parkland Province (6.5 million ha); and Tallgrass Aspen Parklands Province (1.2 million ha) [28].
For the purpose of our study, we focused on fast forest change events, which we define as those that occur within three years or less, which result in immediate changes in spectral indices that are often associated with canopy mortality. Several types of fast forest disturbances occur across the Great Lakes states including land use conversion, timber harvest, wildfire, flooding, extreme weather events, in addition to other natural agents that occur somewhat infrequently (e.g., outbreaks of insect defoliators with immediate impacts on spectral indices). While three years may seem like a long duration for characterizing what we are considering "fast" disturbances, after initial exploration of change events in our study region, we identified the prevalence of multi-year harvest and land conversion activities that we wanted to capture, while the three-year threshold tended to exclude the majority of slower declines events which were not the focus of our study. Land use conversion may take the form of transitioning forest to agriculture, urbanization, or the expansion of mining areas. Timber products are a leading industry in Minnesota, and a range of silvicultural techniques are used that influence harvest intensity, treatment area, duration of activity, and post-harvest site preparation management practices. Harvesting occurs on most lands outside the protected Boundary Waters Canoe Wilderness Area and the Voyagers National Park, with the majority occurring on combined public ownerships compared to private lands. Harvest levels increased substantially from 1970-1995, peaking from 1995-2005, and then decreased to a stable level of approximately 2.8-3.0 million cords annually [29], which represents less than 1% of timberland in the state. Wildfires are most common in the northern part of the state, with prescribed burn events more geographically generalized. Large wind blowdown events are also common in the northern Great Lakes Basin region along with other extreme weather events such as less frequent tornados [1]. Flooding may lead to canopy mortality as a result of prolonged inundation of forested wetland areas, impacts of roadways on drainage, and beaver activity. extreme weather events such as less frequent tornados [1]. Flooding may lead to canopy mortality as a result of prolonged inundation of forested wetland areas, impacts of roadways on drainage, and beaver activity.

Landsat Time Series
We compiled a Landsat time series data set from 1972-2019 to facilitate identification of disturbance events and to aid in the classification of the agent of change. After identifying regional peak growing season characteristics through a sample of annual Normalized Differential Vegetation Index (NDVI) curves, we identified and downloaded Landsat images acquired throughout the growing season from the USGS Earth Explorer website (http://earthexplorer.usgs.gov/) for the years of 1972-2015. Imagery for the years of 2015-2019 were processed within Google Earth Engine (GEE) [30] as a part of later updating efforts described below to develop an efficient work flow for periodic updates of our time series stacks, leveraging recent advances in GEE algorithms and cloud computing. With the high occurrence of clouds in the northern part of the state during the growing season, we attempted to balance the representation of a peak growing season date range with imagery coverage in determining our composite date range, resulting in the time window of July 1st-September 9th. We compiled imagery from all Landsat sensors (MSS, TM, ETM+, and OLI) from 1972-2015 across the entire state of Minnesota, encompassing 28 Landsat scenes with a total of 6041 images included in the initial time series stack. To create annual cloud-free composites of spectral indices that were comparable across all sensors, we utilized LandsatLinkr (LLR) [21] in the R programming environment [31], for all pre-processing and cross-sensor harmonization. LLR automates the necessary steps to create annual composites of the spectral tasseled cap indices (brightness, greenness, and wetness) [32] as well as a derivative of greenness and brightness called the tasseled cap angle [33], harmonized throughout the entire span of the Landsat archive, and composited across multiple Landsat scenes. Incorporating multiple images within the annual compositing time window, including images from overlapping scenes where available, we employed the median compositing method for selecting spectral index values for each pixel across the time series. See Vogeler et al. (2018) for additional details on image processing and LLR steps utilized in this study [11].

Landsat Time Series
We compiled a Landsat time series data set from 1972-2019 to facilitate identification of disturbance events and to aid in the classification of the agent of change. After identifying regional peak growing season characteristics through a sample of annual Normalized Differential Vegetation Index (NDVI) curves, we identified and downloaded Landsat images acquired throughout the growing season from the USGS Earth Explorer website (http://earthexplorer.usgs.gov/) for the years of 1972-2015. Imagery for the years of 2015-2019 were processed within Google Earth Engine (GEE) [30] as a part of later updating efforts described below to develop an efficient work flow for periodic updates of our time series stacks, leveraging recent advances in GEE algorithms and cloud computing. With the high occurrence of clouds in the northern part of the state during the growing season, we attempted to balance the representation of a peak growing season date range with imagery coverage in determining our composite date range, resulting in the time window of 1 July-9 September. We compiled imagery from all Landsat sensors (MSS, TM, ETM+, and OLI) from 1972-2015 across the entire state of Minnesota, encompassing 28 Landsat scenes with a total of 6041 images included in the initial time series stack. To create annual cloud-free composites of spectral indices that were comparable across all sensors, we utilized LandsatLinkr (LLR) [21] in the R programming environment [31], for all pre-processing and cross-sensor harmonization. LLR automates the necessary steps to create annual composites of the spectral tasseled cap indices (brightness, greenness, and wetness) [32] as well as a derivative of greenness and brightness called the tasseled cap angle [33], harmonized throughout the entire span of the Landsat archive, and composited across multiple Landsat scenes. Incorporating multiple images within the annual compositing time window, including images from overlapping scenes where available, we employed the median compositing method for selecting spectral index values for each  [11].
Following the implementation of LLR to produce harmonized spectral indices across the different Landsat sensors, we employed the IDL-based LandTrendr trend fitting algorithm to smooth year to year spectral noise remaining in the dataset, and to identify spectral trends through time [12]. Fitting trend vertices through Landsat time series data sets is a common approach used to identify areas of change, including both fast events (e.g., fire and harvest) and slower declines (e.g., stress and insects), as well as to characterize recovery trajectories. We employed LandTrendr to fit trends through the composite stacks of each of the tasseled cap indices independently, although the trends for tasseled cap wetness (TCW) were ultimately utilized to identify the timing and change magnitude of the most recent fast disturbance event occurring in a given area. Of the tasseled cap indices, TCW is the most sensitive to changes in forest canopy structure [34]. We utilized Google Earth Engine (GEE) to update the time series through 2019, using comparable pre-processing methods and composite date ranges as those utilized in the processing of the 1972-2015 data base, with the exception of fitting LandTrendr trends to the reduced time frame of 1984-2019 to leverage the available Landsat imagery and LandTrendr algorithm [35] on GEE. The fitted tasseled cap indices for 2015-2019 were then exported from GEE and compiled with the existing 1972-2014 fitted indices time series for the identification of disturbance events and classification of agents. Previous studies have noted increased instability of the first year or two of a Landsat time series data base as these early years lack the benefits of previous years for stabilizing context; similarly, the final year of the time series lacks the context of post-disturbance years [12], which may inflate disturbance areas and cause increased errors in agent classifications. Therefore, we allowed the time series two years to stabilize, and removed both of the ending years in the two time series stacks, 2015 from the initial time series, which was replaced by GEE-based indices and 2019 from the second time series, resulting in a final stabilized time series stack from 1974-2018.
We identified the most recent change events from the time series across the landscape, noting the first year post-change as the year associated with each fast event. Since the focus of this study was to identify patch to stand altering events, where we consider patches as within stand events resulting in large openings in the canopy, we followed the methods of Kennedy et al. (2015) and set a minimum mapping unit of 11 Landsat pixels (slightly more than 1 ha) to determine if a patch was retained in the final change map [15]. The minimum of 11 pixels also serves as a homogeneity check to focus on canopy opening events while excluding smaller spectral change patches, which have a greater probability of being related to spatial misregistration issues.
We used a forest mask to remove non-forest lands to limit change patches to just those occurring within forest. There is potential bias in using a current forest mask to remove non-forest areas in a study that spans four decades because areas that were historically forested but recently converted to non-forest would be missed. To minimize this bias, we created an "if-ever-forest-mask" by stacking annual forest masks from 1973-2018 created by Vogeler et al. (2018) [11] and setting a threshold of three years within a forest class at any point throughout the time series to be considered forest during the study period. While this method reduces the bias of removing stands which may have been converted out of a forest class during the study period, it also represents a liberal classification of "forest" for the current time period and may result in some non-forest change events to be retained in the final product.

Change Agent Attribution
Following the identification of the most recent fast change events across the study area, we extracted information about each of the change areas to create a classification model to predict the agent of change for each independent event. Classification predictor variables included: Landsat fitted spectral indices and magnitudes of change, pre-and post-change fitted spectral and canopy cover values, patch area and shape characteristics, topographic indices extracted from Shuttle Radar Topography Mission (SRTM), information about the management protection status as classified by the Gap Analysis Project (GAP) [36] within their Protected Area Data base-US (PAD-US), and key land  Table 1) [37]. Canopy cover values were extracted from an annual state-wide data base, also created using Landsat time series methodologies [11]. In addition to basic elevation and slope, we calculated topographic indices representing local hillslope positions and a topographic landform using the Relief Analysis add-on tool [38] for ArcGIS 10.3. Utilizing the GAP PAD-US product [36], we assigned each change patch with a management protection level from 1 to 4 to include as a categorical predictor, with 1 representing the highest level of resource extraction protection such as those found in national parks and wilderness areas, and 4 indicative of areas with no or unknown management protection. As anthropogenic disturbances should be minimal in highly protected areas, our hope was that the inclusion of this categorical variable may help account for these differences in likelihoods of disturbance agents within the levels of land protection.
We compiled a training database of 2747 patches with known agents of change using a combination of historic aerial photos, existing spatial databases (e.g., Monitoring Trends in Burn Severity) [39], and expert knowledge compiled during personal communications with local land managers. Sources of historical aerial imagery included Google Earth Pro for 1990-present, and the database available through the Minnesota Historical Aerial Photographs Online website for pre-1990 aerial photos [40]. All photo interpretations were completed by one trained interpreter familiar with the study area to reduce biases between multiple interpreters. All training patches were classified using a 6-class designation: land conversion out of a forest class (CONVERISON), forest fire (FIRE), flooding related spectral and/or canopy changes (FLOOD), timber harvest (HARVEST), wind and weather-related events (WIND), and other rare or unknown natural events (OTHER). During the identification of training patches, we noted that some identified change events contained standing water related to flood events, which may or may not have led to actual canopy mortality. We included both flooding resulting in spectral changes, as well as events that ended in canopy mortality, within our FLOOD classification, as spatial information about flooding events, whether they result in immediate canopy mortality or not, can be of interest to managers.
We created a random forest classification model [41] using the randomForest R package [42] for the 6-class agent of change designation. After compiling all predictor variables for training patches, we conducted a variable reduction and model selection procedure employing the rfUtilities R package [43]. We then applied the final 6-class model to all unknown patches to predict the agent of change for all most recent fast change events across the study area.

Change Detection Validation
We followed the methods presented by Olofsson et al. (2014) [44] and Stehman (2014) [45] for the validation of classification maps to quantify area-weighted accuracy measures and confidence intervals for areas of change vs. no-change as well as for each of the change agent classes. We selected a random set of 725 validation map pixels stratified by the six agent classes as well as a seventh class representing areas of forest classified as no-change. We used a combination of the data sources employed in training interpretations (e.g., historic photos) along with plots of the unfitted spectral index trends to aid in the assessment of each of the validation pixels. We assigned a designation of change/no-change and the observed agent of change, which we then compared to the predicted classification to create validation confusion matrices, calculate class and overall accuracies with confidence intervals, and produce an area-weighted coefficient for the estimation of all areas presented in our results.

Spatial and Temporal Summaries
We summarized spatial and temporal trends within forest ownership boundaries as well as within the four ecological regions of Minnesota: eastern broadleaf forest, Laurentian mixed forest, prairie parkland, and tallgrass aspen parkland. We compared the distribution of the different agents of change, magnitudes of agents, as well as time since disturbance classes, across ownerships to identify any spatial or temporal patterns of disturbances within ownerships. Time since last disturbance, along with the magnitude of change, is often related to the succession stage of a stand; thus, in addition to assessing trends in magnitudes of changes by agents, we also utilized the year of onset for each change event to create the following time since disturbance classes with the assumption that they represent different stages of forest succession and associated structures: 0-9 years since disturbance, 10-19 years since disturbance, 20-29 years since disturbance, 30-44 years since disturbance. Leveraging the previously published annual Minnesota canopy cover maps created using the same base Landsat time series spectral data stacks as those used in our study to identify disturbance patterns [11,46], we translated the magnitude of change from each disturbance event into change in percent canopy cover by subtracting the mean pre-disturbance canopy cover from the mean post-disturbance canopy cover, to provide a management-relevant relative measure of magnitudes across agents and ownerships.

Results
The Landsat time series methodology utilized in this study was successful in identifying the most recent fast disturbance events across Minnesota from 1974-2018 with a change vs. no-change  Table 2). The change detection algorithm exhibited higher omission (14.9% ± 9.3%) than commission errors (1.6% ± 1.9%) for the identification of change patches. Omission errors within the validation efforts were often related to the size and intensity of the change event. Change patches close to the minimum size threshold were less likely to be identified, along with patchy natural events where the disturbance impacts were more spread out and barely constituted a contiguous 11 pixels of change. Similarly, lower intensity events also contributed to omission errors, including sivilcultural-thinning activities, which resulted in only minor reductions in canopy. We did not observe any significant differences in error structures of the change and no-change classification if we removed MSS-related change years, supporting its inclusion to extend the time series stack for the identification of fast patch to stand impacting change events. Table 2. Error matrix of estimated proportions of area for change detection assessment and accuracy estimates from independent validation of 725 pixels of the Minnesota disturbance map for the period 1974-2018. Accuracy uncertainty was calculated as 95% confidence intervals and is reported in parenthesis. We further classified disturbance patches by their agents of change (Table 3) with an overall accuracy of 96.5% ± 1.9% (Table 4). Individual disturbance class user and producer errors varied, but all class accuracies were above 78% (Table 4), with HARVEST and CONVERSION exhibiting the highest user accuracies (95.6% ± 3.8% and 92.9% ± 6.2%, respectively) and FIRE and CONVERSION with the highest producer accuracies (98.8% ± 2.4% and 92.6% ± 8.9%, respectively). The lowest user and producer accuracies were observed for the WIND (80.0% ± 9.0%) and HARVEST classes (78.6% ± 10.9%), respectively. Detailed depictions of example patterns of disturbance ages and the classified agents are presented in Figure 2. One of the most prevalent misclassifications were lower intensity harvest patches, such as those that were thinned through strip harvests or localized selective harvests, incorrectly predicted as wind events. There were also some cases where harvests were predicted as flooding and vice versa, which may be due to simultaneous disturbance events (i.e., increase in minor flooding or other moisture related spectral changes immediately following harvest activities) [47] or were a byproduct of post-harvest site-preparation techniques.   During the time period of the study, 17% of forested lands experienced a stand-altering event, although this proportion was calculated using our liberal if-ever-forest-mask, which may bias this estimate through the overestimation of current forested areas; thus, underestimation of proportion of current forested area that experienced a fast change event. Across all regions, harvest was the most prevalent stand replacing causal agent (80% of total disturbed area) followed by wind and flooding, each representing approximately 7% of total disturbance area. Patterns of change agents across the study region show variations in the dominant disturbance events for major stand alterations within the four ecological regions of Minnesota (Figure 3). The Laurentian Mixed Forest ecoregion, which has the largest proportion of public lands including the Chippewa and Superior National Forests, Boundary Waters Canoe Wilderness Area, and Voyagers National Park, had the most variation in agents of change with the highest proportion of area impacted by wind and fire events compared to the other regions ( Figure 3). The Laurentian Mixed Forest ecoregion also had the highest proportion of forested land and disturbed area in general (95% of the disturbed forest). The Eastern Broadleaf During the time period of the study, 17% of forested lands experienced a stand-altering event, although this proportion was calculated using our liberal if-ever-forest-mask, which may bias this estimate through the overestimation of current forested areas; thus, underestimation of proportion of current forested area that experienced a fast change event. Across all regions, harvest was the most prevalent stand replacing causal agent (80% of total disturbed area) followed by wind and flooding, each representing approximately 7% of total disturbance area. Patterns of change agents across the study region show variations in the dominant disturbance events for major stand alterations within the four ecological regions of Minnesota (Figure 3). The Laurentian Mixed Forest ecoregion, which has the largest proportion of public lands including the Chippewa and Superior National Forests, Boundary Waters Canoe Wilderness Area, and Voyagers National Park, had the most variation in agents of change with the highest proportion of area impacted by wind and fire events compared to the other regions ( Figure 3). The Laurentian Mixed Forest ecoregion also had the highest proportion of forested land and disturbed area in general (95% of the disturbed forest). The Eastern Broadleaf Forest ecoregion, home to the densely populated metropolitan area of the Twin Cities, has the largest proportion of private ownership in the state and the highest proportion of land conversion across the study time period (Figure 3).  Land ownerships exhibited varying patterns of likelihood of disturbance across the different change agent classes. The highest proportion of most of the natural disturbance classes, including fire and wind, occurred on federal lands (86% of all fire activity and 50% of all wind activity), which includes the protected Boundary Waters Canoe Wilderness Area and the Voyagers National Park. While harvest was the most likely stand altering disturbance agent across all ownerships (Table 5), the majority of harvest activity in terms of actual land area occurred on private and state managed lands encompassing 469,198 ha (38% of all harvest activity) and 311,553 ha (25% of all harvest activity), respectively. It should be noted that the private and state ownership classes also represent the two largest forest land holdings in the state (Figure 1). Land conversions were most likely to occur on private lands (84% of all conversion area). Land ownerships exhibited varying patterns of likelihood of disturbance across the different change agent classes. The highest proportion of most of the natural disturbance classes, including fire and wind, occurred on federal lands (86% of all fire activity and 50% of all wind activity), which includes the protected Boundary Waters Canoe Wilderness Area and the Voyagers National Park. While harvest was the most likely stand altering disturbance agent across all ownerships (Table 5), the majority of harvest activity in terms of actual land area occurred on private and state managed lands encompassing 469,198 ha (38% of all harvest activity) and 311,553 ha (25% of all harvest activity), respectively. It should be noted that the private and state ownership classes also represent the two largest forest land holdings in the state (Figure 1). Land conversions were most likely to occur on private lands (84% of all conversion area). Assessments of the relative changes in canopy cover within events of different disturbance agents and ownerships revealed significantly higher changes in cover as a result of harvests and land conversions over the other agents, and a tendency for generally higher impacts of events on canopy cover within the "county" ownership class (Figure 4). When we focused on harvest events specifically, the trend for cover magnitudes were similar as those from combined disturbance events, with "county" lands exhibiting significantly higher changes in canopy cover from harvest events (e.g., more clear-cuts) followed by "other public" and "state" managed lands, although harvests within all ownerships spanned a range of sivilcultural techniques and harvest intensities. We acknowledge that the canopy cover maps represent all functional woody vegetation cover over two meters [11], including large woody shrubs common to Minnesota forests, and therefore is not just a measure of canopy removal from the disturbance event, although we believe that it is appropriate as a relative measure of change across the disturbance events. In general, magnitudes of cover changes appeared low, even within canopy clearing events such as clear-cuts, which could be a product of revealed understory trees, a quick release in growth of tall woody understory shrubs, or overestimation of base cover within the annual cover maps. We also calculated change in canopy cover as the mean pre-disturbance cover subtracted from the mean post-disturbance cover within each individual change event; therefore, the resulting change value is impacted by minimal cover change along edges and in areas of retained trees, as well as smoothing over individual pixels within the change patch that may have experienced significantly higher changes in percent cover.  7.02 0.84 Assessments of the relative changes in canopy cover within events of different disturbance agents and ownerships revealed significantly higher changes in cover as a result of harvests and land conversions over the other agents, and a tendency for generally higher impacts of events on canopy cover within the "county" ownership class (Figure 4). When we focused on harvest events specifically, the trend for cover magnitudes were similar as those from combined disturbance events, with "county" lands exhibiting significantly higher changes in canopy cover from harvest events (e.g., more clear-cuts) followed by "other public" and "state" managed lands, although harvests within all ownerships spanned a range of sivilcultural techniques and harvest intensities. We acknowledge that the canopy cover maps represent all functional woody vegetation cover over two meters [11], including large woody shrubs common to Minnesota forests, and therefore is not just a measure of canopy removal from the disturbance event, although we believe that it is appropriate as a relative measure of change across the disturbance events. In general, magnitudes of cover changes appeared low, even within canopy clearing events such as clear-cuts, which could be a product of revealed understory trees, a quick release in growth of tall woody understory shrubs, or overestimation of base cover within the annual cover maps. We also calculated change in canopy cover as the mean pre-disturbance cover subtracted from the mean post-disturbance cover within each individual change event; therefore, the resulting change value is impacted by minimal cover change along edges and in areas of retained trees, as well as smoothing over individual pixels within the change patch that may have experienced significantly higher changes in percent cover.  Within the time since disturbance classes assessed in this study (≤44 years since disturbance), the majority of forest lands across all land ownerships fall within the two youngest time since disturbance age bins (0-9 and 10-19 years since disturbance; Figure 5). County and state-owned lands exhibited higher proportions of areas within the youngest time since disturbance class, while the second youngest class represented the highest proportion of area within federal, tribal, and private lands ( Figure 5). Due to the small forest area within the "other public" category, we did not further evaluate age class distributions within this ownership. Harvest was the driving agent for all time since disturbance classes, with the largest proportion of harvested area within the second youngest class (10-19 years since disturbance; Figure 6). Fire and conversion were more frequent within the youngest class (0-9 years since disturbance), while wind and the other class exhibited their highest influence within the second youngest time since disturbance class ( Figure 6). Within the time since disturbance classes assessed in this study (≤44 years since disturbance), the majority of forest lands across all land ownerships fall within the two youngest time since disturbance age bins (0-9 and 10-19 years since disturbance; Figure 5). County and state-owned lands exhibited higher proportions of areas within the youngest time since disturbance class, while the second youngest class represented the highest proportion of area within federal, tribal, and private lands ( Figure 5). Due to the small forest area within the "other public" category, we did not further evaluate age class distributions within this ownership. Harvest was the driving agent for all time since disturbance classes, with the largest proportion of harvested area within the second youngest class (10-19 years since disturbance; Figure 6). Fire and conversion were more frequent within the youngest class (0-9 years since disturbance), while wind and the other class exhibited their highest influence within the second youngest time since disturbance class ( Figure 6).

Discussion
As forest management agencies increasingly engage in all lands management, spatial information about forest age, condition, and disturbance histories across the landscape may provide a valuable tool for cross-ownership planning. Currently in the western USA, all lands management

Discussion
As forest management agencies increasingly engage in all lands management, spatial information about forest age, condition, and disturbance histories across the landscape may provide a valuable tool for cross-ownership planning. Currently in the western USA, all lands management efforts are becoming more common for fuel reduction treatments in the mitigation of wildfire risks and/or restoration efforts to promote forest resilience to fire events [7], where forest structural stage and disturbance history information may aid in project planning. On a more global stage, studies have begun to explore the utility of cross-ownership assessments for greenhouse accounting frameworks requiring spatially explicit information about forest change, specifically anthropogenic management activities and land conversions, to serve as important accounting and monitoring inputs [48]. More specific to Minnesota, a recent report released by the Minnesota Forest Resources Council identified a key research priority for multi-ownership management related to determining relationships between the spatial arrangement of forest conditions (age class structure, cover types, patch size distribution) and forest resources such as fiber availability, water quality, and wildlife habitat, among others [49]. The approach we have presented here directly addresses information needs of all lands management efforts, and our results show promise for statewide assessments of forest structure and change and its influence on forest resources. With contiguous coverage across the state, independent of land ownership, spatial and temporal information about disturbance patterns presented in our study will help fill information gaps for support of collaborative management and restoration efforts across ownership boundaries within Minnesota, as well as providing a methodological framework, which could be applied to produce similar spatial products in other regions.
In our study, we harmonized imagery across all Landsat sensors, including the early MSS years, which are rarely incorporated into recent Landsat time series change detection analyses. This approach allowed us to identify and classify an extra decade of events compared to many of the other change detection studies [2,15,16,50], resulting in the characterization of forest disturbances across the state of Minnesota for a time period spanning forty-four years. While we did not detect any variations in accuracy structures across the time periods of the different Landsat sensors within our study, we acknowledge that there still may be some sources of limitations and biases in the incorporation of the MSS era not reflected in the validation accuracies alone. While the pre-processing and harmonization procedures facilitated by LLR may have helped alleviate difficulties faced by the coarser resolutions, reduced spectral information, and increased geo-registration errors of the MSS data (see Vogeler et al. 2018 for details on MSS processing steps) [11], there was still the need for extensive quality assessments and periodic manual intervention throughout the processing steps to ensure the creation of comparable products to that of the later sensors. These interventions included the identification and removal of individual MSS images, which caused artifacts in initial spectral composites, as well as manual geo-registration efforts for an additional 96 MSS images (L1G products) to fill in spatial and temporal gaps due to limited higher quality level 1 products, for which LLR is set up to ingest, for some years and Landsat scenes [11]. In addition, the black and white historic photographs from the MSS era used for validation often had larger time gaps between images and coarser resolutions than more recent aerial photography archives, resulting in lower confidence calls during validation assessments. This means that there is a chance that there were additional small lower intensity disturbances that were missed in the MSS years of our disturbance mapping that we were unable to be detected from the available reference data during map validation, which would result in slightly inflated MSS era related accuracies. Other constraints in using MSS imagery include (1) limited availability of MSS imagery and methodologies in GEE, requiring more local computing resources and processing times, and (2) limited coverage of MSS imagery (as readily useable level 1 products) in many areas outside of the U.S., which restricts international opportunities to construct longer time series with MSS data.
Previous studies evaluating time series change detection from 1984-forward utilizing various trend fitting algorithms exhibited similar accuracies to our overall change detection accuracy of 97.2% ± 1.9% [16,50]. Nguyen et al. (2018) utilized a similar LandTrendr segmentation change detection methodology to us, but relied on trends across the normalized burn ratio (NBR) index as opposed to our tasseled cap wetness index for change detection [16]. They reported an accuracy of 87.5% for primary change events across four Landsat scenes within the state of Victoria, Australia [16]. Focusing on the entirety of the province Saskatchewan, Canada, Hermosilla et al. (2015) utilized a Landsat time series break-point detection approach with an overall accuracy of 92.2% for identifying forest change events [50]. In our study, we found higher omission than commission errors for detecting areas of forest change (14.9% ± 9.3% and 1.6% ± 1.9%, respectively), with many of the omission errors attributed to either small (e.g., close to the 11-pixel minimum mapping unit) or lower intensity disturbances (e.g., thinning harvests with high levels of retained canopy). While an underestimation of disturbed area is assumed due to higher omission than commission errors, further validation of the size of individual mapped events could be conducted by digitizing disturbance patches using high resolution reference imagery, if such data exist. Managers and other users of disturbance mapping products may benefit from more detailed validation analyses related to sizes of, and sensitivities to, individual events; thus future work should continue to explore the potential bias of change detection approaches for underestimation in sizes of individual change events and methodological sensitivities to various levels of disturbance intensities (e.g., a range of sivilcultural techniques).
Attributing the causal agent or pattern of change has become a common theme in Landsat time series change detection studies [2,15,16,50] and provides valuable information for land managers and for additional ecological assessments of disturbance regimes and impacts on ecological services. Despite variations in disturbance agent and temporal regimes across regions and forest types, we were able to predict causal agents with class accuracies above 78%, with the majority above 90%, and an overall 7-class (six disturbance classes and one non-change class) accuracy of 96.5% ± 1.9% across the diverse forest systems of Minnesota. When relying on spectral trends for change detection and classification of agents, there are some potential limitations of the data to differentiate between disturbance classes when they result in similar post change structures. For example, patches that experienced lower intensity silvicultural practices such as thinning and selection harvests were occasionally misclassified as wind or missed in disturbance detection altogether. Further, we chose to focus our identification of disturbance areas on relatively fast moderate-high magnitude change events; thus, our maps do not reflect areas impacted by forest changes related to slower declines, which may be the product of insects and disease, although these evaluations are also possible through Landsat time series approaches [51].
Forest harvest represented significantly higher mean changes in canopy cover compared to the other disturbance agents, with the highest mean cover changes occurring on "county", "other public", and "state" managed forest lands. Land ownership can have large impacts on forest disturbance regimes, including sivilcultural intensities, rotation lengths, fire suppression, and salvage harvesting of natural events. For instance, natural events, including fire and wind, were most likely to occur on federal lands. This may be in part due to the historical disturbance patterns of ecological regions, where the majority of federal lands, as well as historic occurrences of large fire events, occur within the Laurentian Mixed Forest Province of Minnesota. The higher likelihood of the most recent change being a wind or fire event on federal lands may also be attributed to the federally owned Boundary Waters Canoe Wilderness Area and the Voyagers National Park, which are protected from harvest and land conversion activities. While wind and fire activities within this ecoregion of Minnesota may also occur on private, state, and county lands, they also have a higher probability of experiencing salvage logging after the natural disturbance event, which would then be reported as the most recent change event within our disturbance map.
By delineating the agent(s) of forest change, we are able to provide greater insight into the residual forest structures, recovery patterns, and potential impacts on ecosystem services. For example, in some instances fire and wind may result in similar canopy mortality as harvest activities, but that reduction in live biomass is likely retained at some level as standing or downed dead wood as opposed to the reduction in canopy being the result of timber removal during silvicultural practices. Those residual standing or downed dead wood legacies within fire or wind impacted stands may provide valuable habitat resources to a variety of wildlife species [52,53] and contribute to the overall structural diversity of the stand. In addition to providing direct insights into residual structures through the identification of change agents, our methodologies may also aid additional modeling and mapping efforts of more detailed structural patterns and stand boundaries through the fusion of remote sensing data sources [19,20]. Combining the historic Landsat time series products with the structural sensitivity of active remote sensing (e.g., lidar) may assist in characterizing difficult to map forest attributes such as biomass [20] and standing dead wood habitat resources [19], which provides value for carbon accounting [54] and monitoring of wildlife habitat for species of conservation interest [53]. Recovery trajectories may also be impacted by the agent of forest change [2]. A study utilizing Landsat spectral trends to assess short-and long-term recovery patterns following fire and harvest events across Canadian forests found that harvest areas recovered more consistently and quickly than fire impacted forests when normalized for magnitudes of change [2]. Through the identification of the timing, extent, magnitude, and agent of disturbance events, we can provide insight into potential recovery trajectories and aid in further research as to the varying impacts of change agents on residual structures and recovery patterns across different forest systems.

Conclusions
As access to affordable remote sensing imagery and cloud computing resources improve, and the temporal records of field and remote sensing-based data bases continue to lengthen, opportunities are increasing to provide spatial records of forest change and recovery patterns across large extents. The findings we presented here demonstrate an approach to quantify most recent fast disturbances over long time periods to provide spatial and temporal information about forest change events and their associated causal agents. These types of data may facilitate cross-ownership monitoring, planning, and restoration efforts to better address the scale of forest ecosystem processes, and improve the efficiency of efforts to address the complicated nature of stressors that forest currently face. The potential of altered disturbance regimes as well as intensities and extents of natural change events with climate change, require drawing on new technologies, management approaches, and large-scale planning efforts to maintain general forest health and functioning across the landscape and the vital ecosystem services that forests provide. The disturbance map presented in this study is freely available on the Minnesota Geospatial Commons (https://gisdata.mn.gov/dataset/env-fast-forest-disturbances).
The unmatched nature of the Landsat archive for providing comparable forest attribute and change information across more than four decades highlights the value of the totality of the Landsat program to the larger geospatial, ecological research, and forest management communities. Our study demonstrates an approach for the harmonization of all Landsat sensors and the identification of change events across the Landsat archive. However, the incorporation of MSS imagery to extend the time series analyses required extensive time and effort, and spatial applicability of our methods are restricted to geographic regions with adequate temporal and spatial coverage of MSS imagery available as terrain corrected level 1 products, which may be limited in some international locations. Further studies should continue to evaluate the values and potential limitations for applications of the extended Landsat time series data set.