Landsat Time-Series for Estimating Forest Aboveground Biomass and Its Dynamics across Space and Time: A Review

: The free open access data policy instituted for the Landsat archive since 2008 has revolutionised the use of Landsat data for forest monitoring, especially for estimating forest aboveground biomass (AGB). This paper provides a comprehensive review of recent approaches utilising Landsat time-series (LTS) for estimating AGB and its dynamics across space and time. In particular, we focus on reviewing: (1) how LTS has been utilised to improve the estimation of AGB (for both single-date and over time) and (2) recent LTS-based approaches used for estimating AGB and its dynamics across space and time. In contrast to using single-date images, the use of LTS can beneﬁt forest AGB estimation in two broad areas. First, using LTS allows for the ﬁlling of spatial and temporal data gaps in AGB predictions, improving the quality of AGB products and enabling the estimation of AGB across large areas and long time-periods. Second, studies have demonstrated that spectral information extracted from LTS analysis, including forest disturbance and recovery metrics, can signiﬁcantly improve the accuracy of AGB models. Throughout the last decade, many innovative LTS-based approaches for estimating forest AGB dynamics across space and time have been demonstrated. A general trend is that methods have evolved as demonstrated through recent studies, becoming more advanced and robust. However, most of these methods have been developed and tested in areas that are either supported by established forest inventory programs and / or can rely on Lidar data across large forest areas. Further investigations should focus on tropical forest areas where inventory data are often not systematically available and / or out-of-date.


Introduction
Landsat satellites are unique in that they have created the longest continuously acquired, consistent space-based and moderate-resolution data collection since 1972. To date, nine Landsat missions have been developed and operated by the United State Geological Survey (USGS). While the Landsat missions 1-4, launched during 1972-1982, acquired data for several years (5-10 years), Landsat 5 delivered high quality, global data across the Earth's land surface for nearly 29 years (1984-2013). Following the failure of Landsat 6 in 1993, Landsat 7 and Landsat 8 were launched in 1999 and 2013,

1.
How has LTS been utilised to improve the estimation of AGB? 2.
What LTS-based approaches have been demonstrated as useful for estimating AGB and its dynamics across space and time?
We first provide a brief review of preprocessing and forest change detection methods for LTS and then focus on answering the above two questions. The first question is addressed by reviewing studies using LTS for estimating AGB, for both a single-date and multiple-dates. The second question is answered by reviewing only studies using LTS for estimating and characterising AGB and its dynamics over time.

Advanced Preprocessing and Change Detection Methods for LTS
Recent studies have included comprehensive reviews on preprocessing and change detection methods for LTS [13][14][15]. For example, Hansen and Loveland [15] reviewed methods for processing LTS for large area land cover change monitoring. Zhu [13] reviewed frequencies, preprocessing and algorithms for change detection using LTS. In this section, we briefly summarise recent methods used for LTS preprocessing and forest change detection that facilitate estimation of AGB and its dynamics.

Robust Preprocessing Methods
Throughout the last decade, many advanced methods have been developed for preprocessing LTS data in order to provide high quality inputs for forest modelling processes. Preprocessing a LTS often includes the following steps: radiometric, geometric and atmospheric corrections, clouds/shadows detection and masking, inter-sensor harmonization, and image compositing. Landsat Level 1 terrain-corrected images have been often used to form a time-series [13]. Radiometric and geometric corrections have been performed on these products, making them suitable for time-series analysis [16]. During 2016-2018, the Landsat archive was reorganised into a formal tiered data collection structure to effectively support long-term LTS "stacking" and analysis (over 40 years) [17]. The new structure manages all Level 1 products in a consistent archive with known data quality. All Landsat acquisitions are stratified into three categories based on data quality: Tier 1, Tier 2, and Real-Time. Tier 1 includes Level 1 Precision and Terrain (L1TP) corrected images that have the highest geometric and radiometric quality (an image-to-image registration accuracy of ≤12-m Root Mean Square Error (RMSE), USGS [17]). Images in this Tier are considered the most suitable for time-series applications. Unfortunately, a large amount of Landsat 1-5 Multispectral Scanner (MSS) images cannot achieve the Tier 1 geometry specification due to less accurate orbital information, insufficient ground control, and other factors including a lack of MSS observations [18]. As a result, LTS relying on Tier 1 products may contain gaps due to the unavailability of high-quality MSS images across some areas and time periods. To address this issue, some studies improved the initial geo-registration of MSS images in Tier 2 (i.e., L1TP with RMSE > 12 m, Systematic Terrain (L1GT) and Geometric Systematic (L1GS) corrections) using some external procedures, which may include some manual steps [19].
Traditionally, atmospheric correction has been conducted using various approaches such as: relative normalization (adjusting the radiometry of other images based on a reference image) [20], the Dark-Object Subtraction (DOS) [21], and the Satellite Signal in the Solar Spectrum (6S) [22]. However, these methods often require many manual and complex steps which hinder the processing of large amount of LTS [13]. Recently, two algorithms are being implemented for automatic atmospheric correction: the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) [23] for Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+), and the Land Surface Reflectance Code (LaSRC) [24] for Landsat 8 Operational Land Imager (OLI) images. LEDAPS and LaSRC have been used to generate surface reflectance images from Landsat Level 1 products that are currently provided on the USGS archive for public access (identified as Landsat Level 2 products). These surface reflectance products have been increasingly used in recent LTS applications as they are robustly pre-processed [13]. Another benefit of using surface reflectance images is that they are provided with high quality cloud/shadow detection products. Cloud/shadow in LEDAPS and LaSRC products are detected using the CFMask algorithm, originally known as FMask developed by Zhu and Woodcock [25]. This algorithm is widely confirmed to provide a higher accuracy of cloud/shadow detection than other methods. While surface reflectance and cloud/shadow masks for TM, ETM+ and OLI images are consistently created and available on the USGS archive, additional preprocessing steps are required to derive these data from MSS images. MSS surface reflectance images are often calculated using the COST method [26], while cloud/shadow can be identified using the MSS clearview-mask (MSScvm) algorithm [27].
A challenge of LTS analysis is harmonising the different spectral, spatial, and radiometric resolutions of Landsat sensors. Most LTS studies utilise images acquired by TM and ETM+ sensors (1984-present) for their time-series analyses since they have similar spatial and spectral characteristics [19]. Utilising the full temporal record of Landsat data (including MSS and OLI data) is valuable for characterising forest dynamics, but it requires additional processing to normalise and incorporate multi-sensor images into a consistent time-series [12,28]. Several inter-sensor harmonization approaches have been recently demonstrated to expand the utility of the Landsat archive for forest dynamics applications. Vogeler, Braaten [19] presented an automated method, called LandsatLinkr, for calibrating Landsat MSS and OLI images to the spatial and spectral profiles of TM and ETM+ images. LandsatLinkr normalises MSS images to the TM images by modelling TM Tasseled Cap (TC) indices from the four MSS surface reflectance bands using multiple linear regression. Regression models are developed based on the relationship between coincident MSS and TM images simultaneously acquired by Landsat 4-5 between 1982 and1992. They are then applied to normalise all MSS acquired by Landsat 1-3 (back to 1972). A similar approach is deployed in LandsatLinkr for harmonising OLI images to ETM+ images.
Given the limited availability of cloud-free Landsat data in many forested areas, image compositing is an important process for supporting AGB predictions, especially temporal predictions (see Section 3). This process leverages the extensive LTS data to create cloud-free and radiometrically consistent composite images that are spatially contiguous over large areas. Image compositing is also useful for reducing data volume and minimizing atmospheric influences. Since the opening of the Landsat archive in 2008, several pixel-based approaches have been developed for generating composite images from Landsat data. A robust method commonly used in LTS applications is the Best Available Pixels (BAP) compositing, developed by White, Wulder [29]. This method determines the BAP within a compositing window/interval (e.g., annual, season or month) using pixel-scoring functions including pixel-based scoring for sensor, acquisition date, distance to cloud/shadow, and atmospheric opacity [30]. Zhu, Woodcock [31] demonstrated an approach for generating daily clear-sky synthetic images from all available Landsat data by developing time-series models for each pixel and each spectral band. Other approaches such as using maximum, mean or median values of Landsat bands and indices are also used for Landsat image compositing [32][33][34].

Vegetation Change Detection Using LTS
Having knowledge of the spatial and temporal patterns of forest disturbance and recovery is crucial to improving the accuracy of AGB estimates and associated dynamics (Sections 3 and 4). Forest disturbances can be abrupt or subtle events caused by natural hazards such as wildfires, drought, and insects or human activities such as logging and mining [35]. In contrast, forest recovery is an ongoing process that can occur naturally through vegetation succession or artificially through forest management [36]. The ability of LTS to capture changes occurring in forests has been widely demonstrated, with various methods developed to characterise forest disturbance and recovery. In a recent review, Zhu [13] categorized LTS change detection methods into six broad groups based on the algorithms used for detecting changes (e.g., differencing, thresholding, trajectory classification, regression, statistical boundary and temporal segmentation). Of these, thresholding, statistical boundary and temporal segmentation have been the most commonly used. Table 1 lists some popular algorithms based on those three methods. These pixel-based algorithms have proven to be capable of detecting changes across large areas and long time-periods. In addition, they are fully automated, publicly available and well documented algorithms. The VCT is a threshold-based algorithm. It calculates IFZ as a forest probability index for each Landsat pixel in a time-series and determines an IFZ threshold to detect forest disturbance (pixels with an IFZ value smaller than the defined threshold are considered as disturbance) [37]. VCT provides information about abrupt disturbances on an annual scale using an annual LTS. In contrast, CCDC and BFAST Monitor estimate a statistical boundary for the entire time-series and define changes as significant differences from the boundary [25,39]. Both of these algorithms utilize all available Landsat images and are designed to detect both abrupt and gradual forest change. NDVI is often used with BFAST Monitor while CCDC runs on all Landsat spectral bands. LandTrendr is a temporal segmentation algorithm utilising an annual LTS. This algorithm fits the spectral trajectory of each pixel to a series of straight-line segments to capture multiple change events as well as the overall trend of the trajectory over time. LandTrendr is capable of detecting both abrupt and gradual changes at an annual timescale. Although NBR is often used as the change index for LandTrendr, other vegetation indices (such as TC components and NDVI) and spectral bands have also been deployed [41,42]. The robustness of LandTrendr has been demonstrated across various contexts including different geographical regions, forest types and disturbance regimes [13,14]. The algorithm has been recently implemented on the Google Earth Engine platform, offering a consistent and quick approach for mapping forest disturbance and recovery from LTS data [43].
The utility of a change detection technique is dependent on its ability to detect and distinguish real changes (or changes of interest) from apparent changes caused by variations in background signal or noise [42]. Studies have demonstrated that high-impact or stand-clearing disturbances and subsequent recovery can be accurately detected from LST analysis. Whereas, changes caused by medium/low-impact disturbances (such as selective logging and droughts) are more difficult to distinguish and characterise [44]. Therefore, the selection of the most appropriate change detection algorithm for a specific application is challenging and dependent on multiple factors including target changes and the temporal density of LTS used. VTC mainly targets abrupt, high magnitude forest disturbances while CCDC and LandTrendr target a broader range of forest changes (both discrete and gradual changes, and both disturbance and recovery processes). Generally, methods requiring a higher density of LTS (i.e., CCDC and BFAST Monitor) are often computationally expensive and require long-processing times and substantial data storage. In many instances, an annual LTS should be the most suitable temporal scale for mapping forest changes as it allows for the tracking of both abrupt and gradual change without requiring as much processing effort. But regardless of the choice, a benefit of using automate change detection algorithms is that they can robustly produce a suite of spectral change metrics representing trends of forest disturbance and recovery. Change metrics are often utilised in many applications such as attributing disturbance levels and causal agents and improving forest biomass estimates [12,[45][46][47].

How Has LTS Been Utilised to Improve the Estimation of AGB?
Since the opening of the Landsat archive in 2008, LTS data have been increasingly used for estimating forest AGB over large areas, and not only for a static date but also for multiple points in time. The contribution of LTS to improvements in forest AGB estimates has been demonstrated in a wide range of contexts. In Table 2, we summarise the benefits of using LTS for forest AGB estimation into two broad categories: (1) filling spatial and temporal data gaps in AGB predictions and (2) improving the accuracy of AGB modelling. Filling spatial and temporal data gaps in LTS using a pixel-based temporal fitting process Deo, Russell [57]; Deo, Russell [58]

Benefit Specific Improvement/Finding Reference
Improving the accuracy of AGB modelling Using disturbance and recovery metrics derived from LTS can improve the accuracy of AGB predictions in comparison with using single-date images Pflugmacher, Cohen [28]; Pflugmacher, Cohen [12] Quantifying the capacity of LTS-derived change metrics for estimating forest AGB Frazier, Coops [54] Demonstrating the utility of LTS-derived change attributions (e.g., disturbance severity and agent) in modelling forest structure and AGB Zald, Ohmann [59]; Zald, Wulder [48]; Nguyen, Jones [56]; Bolton, White [60] Forest age derived from LTS can improve forest AGB estimates Lefsky [61]; Liu, Peng [62] Identifiers of temporal patterns in spectral trajectories provide improvements in modelling and explaining historical AGB dynamics Gómez, White [63] Seasonal NDVI time series improving AGB estimations compared with a single-date NDVI Zhu and Liu [64] Fitted LTS data improving the accuracy of AGB predictions compared with raw/observed data. Deo, Russell [58] A pixel-based temporal fitting process improves the temporal consistency of AGB predictions.

Filling Spatial and Temporal Data Gaps in AGB Predictions
Forest AGB is commonly estimated using high quality Landsat images (i.e., cloud-free observations) from small areas e.g., [65][66][67][68][69][70][71][72][73][74]. However, when larger areas and long time-periods are desired, it can be challenging to find cloud-free Landsat data, resulting in spatial and temporal gaps in AGB prediction maps. This prohibits creation of wall-to-wall AGB maps over large areas (i.e., across multiple Landsat scenes) and/or multiple points in time. LTS data provide an opportunity to overcome these limitations by computing image composites (Section 2.1). Image compositing methods rely on selecting representative pixels (e.g., via BAP) from a series of images rather than the best available scene [29]. Landsat image composites are increasingly used to support the development of spatially and temporally complete predictions of forest AGB (Table 2). Zald, Wulder [48] indicated that Landsat pixel-based composites are important for improving the quality of AGB predictions over large areas, regardless of the modelling technique and training data used. Matasci, Hermosilla [49] demonstrated the utility of Landsat composites for extending localised concept studies to the national scale by creating spatially explicit estimates of AGB across 552 million ha of Canadian boreal forests. The authors also extended their work to estimate AGB dynamics across three decades (1984-2016) [51]. Despite the significant increase in the number of cloud-free observations, following the launch of Landsat 8 in 2013 (and Landsat 9 expected in 2020), LTS-based image composites will remain essential for forest AGB monitoring across large areas and long-time periods, especially for regions with persistent cloud and snow cover [29,48].
Pixel-based temporal fitting is another process that can be performed on LTS to fill (or further fill) spatial and temporal data gaps [48,58,75]. This process results in a "smoothed" or "stabilized" temporal trajectory for each Landsat pixel, filling data gaps (i.e., no data observations) and removing noisy pixels (due to unscreened cloud/shadow). Pixel-based temporal fitting can be used for filling not only spatial gaps but also temporal gaps resulting from the unavailability of Landsat data at some points in time. For example, Deo, Russell [58] applied a pixel-level curve fitting approach to a LTS of 17 cloud-free images representing 26 years (1986-2011) to obtain fitted images for years where cloud-free images were not available. Pixel-base temporal fitting for a LTS can be conducted using curve fitting tools [58] and temporal regression or segmentation algorithms such as LandTrendr [9,52].

Improving the Accuracy of AGB Modelling
While a pixel-based temporal fitting process can be used for filling data gaps, a more important purpose of this analysis is to minimize temporal noise/variations (due to exogenous factors such as sun angle, phenology and atmosphere) at the pixel-level of LTS [9,51,[57][58][59]76,77]. It is demonstrated that this process enables the consistent implementation of a statistical model over time and improves the accuracy of biomass predictions ( Table 2). The use of calibrated predictor values is the key factor enabling the consistent implementation of a statistical model across both space and time [20]. Advanced developments in image pre-processing (Section 2.1) facilitate the generation of normalized time-series images with minimal data gaps and a consistent radiometric response. This established a foundation for estimating forest AGB measurements through time [51]. However, modelling continuous forest attributes such as AGB can be negatively impacted by residual temporal noise or variations existing in Landsat pixel-level spectral trajectories [28,48,57,58]. Thus, a pixel-based temporal fitting process is necessary to mitigate temporal fluctuations and thus improve the temporal coherence in AGB predictions. Deo, Russell [58] demonstrated that using spectral data (i.e., spectral bands and indices) extracted from a fitted LTS can improve the accuracy of AGB predictions (RMSE = 47.64 Mg·ha −1 ) in comparison with using observed data (RMSE = 54.18 Mg·ha −1 ), as the temporal noise has been removed through a fitting process.
Pixel-based temporal fitting analysis captures spectral changes and trends of forested pixels (due to forest disturbance and recovery processes) over time [37,38,40,42,78]. It is widely confirmed that LTS-derived change metrics characterising the patterns of forest disturbance and recovery can improve the accuracy of AGB estimates ( Table 2). An early example of this application showing the utility of disturbance and recovery metrics derived from an annual LTS analysis  in predicting forest AGB is presented by Pflugmacher, Cohen [28]. Spectral change metrics such as disturbance and recovery onset years, duration, magnitude, and time since disturbance (TSD), pre-and post-change conditions were derived from a temporal segmentation analysis using the LandTrendr algorithm [40]. The study compared LTS-based models with traditional modelling approaches which rely on single-date Landsat imagery and Lidar data. They found that forest disturbance and recovery metrics substantially improved the accuracy of model predicting AGB of live trees (R 2 = 0.80, RMSE = 46.9 Mg·ha −1 ) in comparison with single-date Landsat data (R 2 = 0.58, RMSE = 65.1 Mg·ha −1 ). Moreover, LTS-based models were significantly better in predicting AGB of dead trees (R 2 = 0.73, RMSE = 31.0 Mg·ha −1 ) than Lidar and single-date models (R 2 = 0.21, RMSE = 43.8 Mg·ha −1 ; R 2 = 0.00, RMSE = 47.8 Mg·ha −1 , respectively). The capacity of LTS-derived change metrics for estimating forest attributes, including AGB, was further quantified by Frazier, Coops [54]. In this study, a number of disturbance and recovery metrics (both simple and complex metrics) were calculated from LandTrendr-derived spectral trajectories of TC components (TCB, TCG and TCW). Eleven RF models were then tested based on different groups of predictor variables to investigate the importance of change metrics for estimating forest biomass. The study concluded that simple change metrics (such as change onset year, duration and magnitude, pre-and post-change values) have more predictive power than complex metrics (such as longest recovery, year to year change rates, and last monotonic trends). In addition to spectral change metrics, change attributions derived from LTS such as disturbance type (or causal agent) and severity level have been recently used as predictor variables for forest AGB modelling [9,49,52]. These metrics can contribute to the improvement of prediction accuracy by distinguishing pixel-level spectral conditions. For instance, pixels can experience a similar spectral change magnitude, at the same point in time, but different casual agents (e.g., fire and logging), resulting in different forest structure and biomass [47,59]. The utility of other LTS-based change metrics such as forest age has also been demonstrated in different contexts [61,62]. Liu, Peng [62] found that including LTS-derived forest age significantly improved the accuracy of AGB model when compared to using spectral indices only (R 2 value increased from 0.50 to 0.73). The use of LTS change metrics in AGB estimation has been facilitated by developments in change detection algorithms (Table 1) and innovative cloud-based methods, such as the integration of the LandTrendr algorithm on Google Earth Engine [43].
As an exception, Gómez, White [63] used "dynamic" variables extracted from LTS pixel trajectories instead of spectral change metrics. "Dynamic" variables included transformed spectral trajectories and temporal derivatives associated with those trajectories. These variables were derived as identifiers of temporal patterns in spectral trajectories using two-dimensional (2D) wavelet transforms of an ordered system. The study demonstrated that "dynamic" variables provided some improvements in modelling and explaining historical AGB variability. However, the extraction of such "dynamic" variables is dependent on the availability of re-measured AGB data.
Though most applications use an annual or near-annual LTS, some studies have demonstrated the potential of using a higher temporal density LTS to improve the estimation of forest AGB [64,79]. Zhu and Liu [64] explored the use of seasonal LTS and found that seasonal NDVI time-series can improve the accuracy of AGB estimations compared with a single-date NDVI (R 2 = 0.54 vs. R 2 = 0.42). The study also indicated that NDVI in the fall season (September and October) has stronger correlation with AGB (R = 0.64) than in the peak season (June, R = 0.46). Nguyen, Jung [79] used LTS to estimate seasonal AGB for mixed coniferous forests in South Korea. Results from this study indicated that using Landsat images in summer (August) can provide more accurate estimation of forest AGB than using images from other seasons, although seasonal accuracies were relatively stable. The utility of using high temporal LTS (i.e., seasonal or denser scales) should be further investigated across different contexts.

What LTS-Based Approaches Have Been Demonstrated for Estimating AGB and Its Dynamics across Space and Time?
In this section, we review recent studies using LTS data as inputs for estimating forest AGB and its dynamics across space and time (Table 3). Of note, all studies investigated boreal and temperate forests in North America (US and Canada) and Europe, none of them was conducted in tropical forests. Study areas ranged from regional scales (several Landsat scenes) to continental scales such as over 650 million ha of Canadian forests [51]. Commonly investigated periods ranged from 20 to 40 years. Most of the studies used Landsat TM/ETM+ data to form a LTS stack (starting in the 1980s) while few of them utilized MSS data to extend the time-series further back to 1972 [12,62]. Studies often developed LTS stacks at an annual timescale using cloud-free observations or image composites. However, multi-year epochal or seasonal-scale approaches were also used in some applications. In addition, all studies estimated historical AGB, except Ma, Domke [80] who made AGB future projections.    As shown in Table 3, many LTS-based approaches for estimating and characterising forest AGB dynamics have been carried out in different contexts and they have become more advanced and robust over time. A conceptual diagram for estimating AGB dynamics from LTS is shown in Figure 1. LTS data are processed to extract explanatory variables (e.g., spectral indices and change metrics). These are then used in a modelling framework as predictor variables along with reference data (derived from forest inventory or Lidar-based plots). The implementation of this step results in predictions of AGB across space and time. Patterns of AGB dynamics are often characterized by integrating AGB changes with forest disturbance and recovery history, which are normally derived from LTS. Accuracy of predictions is often assessed using independent field plots or Lidar-based AGB data. A few studies, however, followed other approaches (e.g., Ma, Domke [80] and Zhang, Lu [84]) to directly predict AGB dynamics from LTS (Table 3).
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 25 As shown in Table 3, many LTS-based approaches for estimating and characterising forest AGB dynamics have been carried out in different contexts and they have become more advanced and robust over time. A conceptual diagram for estimating AGB dynamics from LTS is shown in Figure  1. LTS data are processed to extract explanatory variables (e.g., spectral indices and change metrics). These are then used in a modelling framework as predictor variables along with reference data (derived from forest inventory or Lidar-based plots). The implementation of this step results in predictions of AGB across space and time. Patterns of AGB dynamics are often characterized by integrating AGB changes with forest disturbance and recovery history, which are normally derived from LTS. Accuracy of predictions is often assessed using independent field plots or Lidar-based AGB data. A few studies, however, followed other approaches (e.g., Ma, Domke [80] and Zhang, Lu [84]) to directly predict AGB dynamics from LTS (Table 3).

Explanatory Data
Explanatory data used for multi-temporal AGB modelling generally include three groups of variables: (1) Landsat image-based spectral variables (i.e., Landsat surface reflectance bands, vegetation indices, band ratios, and texture metrics), (2) LTS-derived spectral change/trend metrics and (3) other ancillary data such as topographic and climatic layers (Table 3).
Image-based spectral variables are key predictors for estimating forest AGB from Landsat data. These tend to incorporate the near-infrared (NIR) and shortwave infrared (SWIR) wavelengths such as TM/ETM+ band 5, OLI band 6, NBR, NDVI and EVI, given they often have strong relationships with observed AGB values [10,85]. The TC transformations, including TCB (overall reflectance), TCG (contrast between NIR and visible reflectance), and TCW (contrast between NIR and SWIR) [86][87][88][89], are also used due to their abilities in revealing key forest attributes, including AGB [28,54,85,90]. In addition to indices directly derived from Landsat bands, the use of LTS has facilitated the introduction of some pseudo/integrated spectral indices that are useful for AGB modelling, such as TCA, TCD, DI and IFZ (see Appendix A for descriptions). These indices are often used with an assumption that using an integrated index may provide more useful information than using a group of individual indices. Healey, Cohen [91] showed that using DI alone provided better results in a land

Explanatory Data
Explanatory data used for multi-temporal AGB modelling generally include three groups of variables: (1) Landsat image-based spectral variables (i.e., Landsat surface reflectance bands, vegetation indices, band ratios, and texture metrics), (2) LTS-derived spectral change/trend metrics and (3) other ancillary data such as topographic and climatic layers (Table 3).
Image-based spectral variables are key predictors for estimating forest AGB from Landsat data. These tend to incorporate the near-infrared (NIR) and shortwave infrared (SWIR) wavelengths such as TM/ETM+ band 5, OLI band 6, NBR, NDVI and EVI, given they often have strong relationships with observed AGB values [10,85]. The TC transformations, including TCB (overall reflectance), TCG (contrast between NIR and visible reflectance), and TCW (contrast between NIR and SWIR) [86][87][88][89], are also used due to their abilities in revealing key forest attributes, including AGB [28,54,85,90]. In addition to indices directly derived from Landsat bands, the use of LTS has facilitated the introduction of some pseudo/integrated spectral indices that are useful for AGB modelling, such as TCA, TCD, DI and IFZ (see Appendix A for descriptions). These indices are often used with an assumption that using an integrated index may provide more useful information than using a group of individual indices. Healey, Cohen [91] showed that using DI alone provided better results in a land cover change analysis than using a group of TC components or TM reflectance bands. Since first introduced by Powell, Cohen [81], TCA, a combination of TCB and TCG, has been utilised as an important predictor for estimating forest AGB dynamics in a wide range of research contexts (Table 3). Though IFZ (representing the likelihood of a pixel being forested) was initially used for forest change detection [37], it has been used for AGB modelling in some recent studies [57,58,80]. According to Deo, Russell [58], this metric provides a normalized predictor that can significantly reduce the spatial and temporal variability of spectral signatures caused by atmospheric conditions and sensor issues, improving the accuracy of AGB predictions. However, the utility of this index in estimating forest AGB has only been demonstrated in deciduous and coniferous forests (North America) and thus needs to be further investigated in other forested areas.
In the context of estimating forest AGB over time, the selection of spectral indices is not only dependent on their relationship with observed AGB values but also their ability to capture the forest's history of disturbance and recovery. Among reviewed studies (Table 3), spectral indices such as NBR, NDVI, TCW and TCA are the most popular given their sensitivity to changes occurring in forests. These indices are often analysed through a temporal fitting process to extract forest disturbance and recovery metrics that contribute to the improvement of AGB predictions (Section 3). Another consideration when selecting spectral indices for estimating forest AGB over time is which sensors are being used. Studies utilising MSS data often relied on simple band ratios [62] or modelled TC indices [12]. Some spectral indices that are useful for AGB estimation and forest change detection such as NBR cannot be calculated from MSS images due to the lack of SWIR bands. This suggests that the inclusion of MSS data in LTS should be specifically considered for each case study. Spectral indices that are the most important for AGB estimation, and also the most sensitive to forest change, should be determined before processing MSS data. If selected indices are not available from MSS data, then the inclusion of MSS images can negatively impact temporal predictions of forest AGB.
Forest AGB estimates can be enhanced by incorporating spectral change metrics representing disturbance and recovery trends prior to the prediction date ( Table 2, Section 3). However, the use of disturbance and recovery metrics for multiple-date AGB predictions is different from a single-date estimation. A single-date estimation can utilise temporal trends extracted from the full length of LTS. In contrast, when estimating AGB over multiple dates, the length of LTS diminishes as the prediction date moves further back in time and consequently, change metrics contain less information on disturbance and recovery which may negatively impact prediction accuracy [12]. Pflugmacher, Cohen [12] found that model accuracy gradually reduced when decreasing the length of LTS. Furthermore, they noticed that the appropriate length of LTS needed to derive meaningful change metrics is dependent on the disturbance frequency and intensity of a given study area (in their case between 10-20 years). Gómez, White [63] demonstrated that 15-25 years is sufficient for identifying significant temporal patterns in spectral trajectories to support AGB estimates in pine forests, Spain. This suggests that the balance between the prediction accuracy and the required temporal density will depend on the specific application. If only periodic dates are required for making predictions (e.g, more than 10 years), spectral disturbance and recovery metrics should be included to optimize the accuracy of biomass predictions. In contrast, if predictions are required at higher temporal densities (e.g., annual), single-date predictor variables should be used to maintain the temporal accuracy. Recent studies have incorporated some change attributions, such as disturbance causal agents and TSD, in annual modelling of forest biomass [9,51,52,62]. Though the predictive power of these metrics is normally not significant, they can benefit AGB modelling by distinguishing pixels with similar spectral information [48].
In addition to LTS-based variables, ancillary data such as topographic and climatic information are often included in AGB estimating across large areas ( Table 3). The importance of these variables for predicting AGB has been widely highlighted given topographic and climatic conditions are often diverse over large areas, resulting in different forest ecosystems.

Training Data
For predicting AGB across space and time, LTS data need to be combined with plot training data. Studies often use two main sources of AGB training data: forest inventory and Lidar-based plots. Observed AGB estimates are often directly calculated from forest inventory plots using species-based allometric relationships between tree metrics and their biomass. On the other hand, Lidar-based AGB plots are predicted through a modelling process that combines Lidar with forest inventory data.
Forest inventory is the most common source of AGB reference data. A significant benefit of using forest inventory plots is that they provide the most accurate observations of AGB. Unfortunately, the number of available field plots is often limited across large and remote areas. When developing a LTS-based model, inventory plot data collected from multiple points in time and data sources can be incorporated, increasing the sample density for AGB modelling. For example, Kennedy, Ohmann [9] integrated all available plots from several inventories surveyed during 1991-2011 into a database of 4318 plots with 8454 measurements. The authors then developed a monolithic model to estimate annual AGB by combining the plot database with LTS data (extracting spectral values from image dates coinciding with plot measurement dates). In another study, Boisvenue, Smiley [50] used 1381 field plots across Saskatchewan, Canada (originating from 1949) and LTS data to estimate AGB dynamics from 1984 to 2012. It is important to note that uncertainties arising from ingesting different data sources need to be carefully considered given the inconsistency in inventory methods including plot sizes and survey techniques [9,50]. The use of forest inventory data is also dependent on modelling methods used for estimating AGB dynamics. Some approaches rely on data from remeasured plots which are often not available in many forested areas [50,80,84].
Airborne Lidar data provide the most accurate remote sensing-based predictions of forest AGB, but they are often not available wall-to-wall and over time given high acquisition costs and computational requirements. Alternatively, recent studies often employ airborne Lidar data as a sampling tool to improve the efficiency of forest inventory data. Sample-based Lidar plots can significantly increase the number of accurate samples of AGB estimates in areas without sufficient field inventory data [92], reducing the variability in plot-level AGB values. Lidar-based plots can be used as an alternative for training models based on LTS data, which are temporally and spatially complete [48,49].
The utility of Lidar-based plots in estimating forest AGB dynamics has been demonstrated in several studies (Table 3). Pflugmacher, Cohen [12] used airborne Lidar to sample a wide range of forest disturbance histories, as an alternative of 51 field plots within their study area. Lidar-based AGB plots were then combined with LTS data for large-area modelling of annual forest AGB over 38 years (1972-2009). Matasci, Hermosilla [51] extracted nearly 85,000 Lidar plots from 34 survey transects across Canadian forest ecosystems and integrated these plots with LTS to produce annual AGB estimates and related dynamics from 1984 to 2016 at the national scale. However, Lidar-based AGB predictions contain uncertainties that should be considered when describing the uncertainties in Landsat-based predictions [12,48]. Although Lidar data often show a strong relationship with observed AGB values, the consequences of introducing additional uncertainty into Landsat-based AGB maps need to be further investigated. Furthermore, Lidar-based plots are often spatially constrained by Lidar transects, which may not cover the full range of the forest population (e.g., transects may not be systematically distributed across different forest types). This will add more uncertainty into AGB predictions, especially in the case of a kNN imputation model, with k = 1, being used (i.e., predicted values cannot exceed the range of sample dataset) [48]. Another limitation is that Lidar-based plots cannot provide species composition information due to the lack of tree measurements; therefore, they cannot be used when deriving species specific AGB estimates.

Statistical Modelling Technique
Various techniques have been deployed for modelling AGB across space and time. These can be grouped into three broad model types: parametric regression, nonparametric regression and nonparametric imputation. Parametric regression is a traditional method for estimating AGB based on remote sensing, which finds the linear relationship between observed AGB and predictor variables. Few studies opt for this method (Table 3), given it is known to overestimate and underestimate low and high AGB values, respectively. Furthermore, Powell, Cohen [81] and Zhang, Lu [84] compared different empirical approaches for modelling AGB dynamics from LTS and both found that regression models, such as Linear Regression and Reduced Major Axis, provide less accurate results than non-parametric methods such as Random Forest (RF) and Gradient Boost Regression Tree.
Non-parametric approaches are increasingly used to derive forest structural attributes from remote sensing data as they do not have restrictive assumptions on the data analysed. Two nonparametric regression approaches commonly used for modelling AGB dynamics are regression trees and kNN [93] ( Table 3). Regression tree methods are widely used for forest applications due to their tremendous analytical and operational flexibility. RF is a regression tree algorithm widely used for forest modelling. It is an ensemble learning algorithm that develops a number of small regression trees that vote on predictions. Each tree is constructed from a random subsample from the original dataset and nodes on trees are split using a random subset of predictor variables [94]. RF is capable of accounting for complex and non-linear relationships among variables, and robust in preventing over-fitting. An advantage of using RF is that it can handle a large number of variables/metrics extracted from LTS data. In addition to RF, other regression tree algorithms such as Decision tree (known as CART) and Gradient Boost Regression Tree have also been used for modelling AGB dynamics [63,84]. kNN regression predicts the AGB value for a target sample (pixel) as the mean of the k most similar reference samples, defined using a distance metric. Although it is a commonly used modelling method for forest applications [95], few studies have used kNN regression for estimating AGB over time based on LTS [79]. Nevertheless, Wilson, Knight [55] demonstrated that kNN regression and RF produced comparable accuracies when modelling forest AGB from LTS data.
Non-parametric imputation approaches have recently received considerable attention for forest inventory mapping applications. These multivariate modelling approaches can simultaneously map large assemblages of inventory attributes, therefore retaining the variance structure to that of the observations [96]. For predicting forest AGB dynamics from LTS, the kNN imputation has been frequently employed (Table 3). This method imputes (or shares) observed values of one or multiple response variables to target samples/pixels [8]. Kennedy, Ohmann [9] demonstrated that the kNN imputation approach performed well to overcome the saturation of Landsat signals known to occur in forests with high biomass. Unlike regression, where a predicted value is often different from any observed values, an imputed value will be the most similar observed value when the single nearest neighbour is used (k = 1). When k > 1, an imputed value will be the mean of observed values associated to the k nearest neighbours. While using a k > 1 can improve the imputation accuracy and produce more smoothed predictions, the use of k = 1 is preferable to maintain the variance structure of predictions that are similar to the observations [76]. Most studies deploying kNN for AGB dynamic mapping used the single nearest neighbour as they preferred to utilize the full range of AGB values in training data [9,51,52]. Deo, Russell [58] suggested that a k value of 3 to 5 is the most suitable in order to improve model accuracy while retaining the variance structure of predictions closer to observations. Nguyen, Jones [52] achieved this using k = 1 when processing the kNN model though a bootstrapping analysis with multiple repetitions.
In kNN imputation, the nearest neighbours are found based on a distance metric derived using different techniques such as Gradient Nearest Neighbour (GNN) [68] and RF [97]. GNN is a weighted Euclidean distance-based technique which is a combination of the canonical correspondence analysis and k = 1. In contrast, RF is a machine learning-based technique in which the distance metric is calculated based on a proximity matrix [98]. Several studies have shown the utility of GNN in modelling AGB dynamics [9,81]. In a recent study, Kennedy, Ohmann [9] integrated LandTrendr and GNN-based imputation in a robust empirical forest biomass monitoring system. RF-based kNN imputation has been increasingly used in various contexts for estimating AGB dynamics based on LTS data. For example, Matasci, Hermosilla [51] and Nguyen, Jones [52] used this method for large area predictions of AGB dynamics through 30-year periods in Canada and southeast Australia, respectively. Studies comparing different kNN distance techniques used for imputing forest attributes often indicated the outperformance of RF [8,56,96].
When using kNN imputation, AGB can be predicted using either a direct or indirect imputation approach. In the direct approach, AGB is imputed by directly linking observed AGB (as a response variable) with Landsat-derived predictors through a kNN model [58,81]. However, recent studies often used an indirect imputation approach for estimating AGB. In this method, a kNN model is trained by forest structural variables extracted from inventory plots [9,52] or Lidar metrics extracted from Lidar plots [51]. AGB is not included in the model and is attached as an ancillary variable to the imputation process for target samples. In other words, the nearest neighbours are found according to the relationship between structural variables and LTS-derived predictors and then indirectly transferred to AGB. Response variables included in the model should appropriately capture forest conditions and can be explained by Landsat spectral data. The accuracy of AGB predictions can be enhanced as it is imputed based on another variable (or group of variables) that has higher correlation with predictors. Nguyen, Jones [56] found that AGB can be better predicted using an indirect imputation approach combining a group of basal area and stem density variables with LTS. In addition, other forest attributes, along with AGB, can be mapped directly without new model developments. Kennedy, Ohmann [9] combined field plot measurements, including basal area by species and tree size classes, with LTS-derived predictors through a GNN-based imputation model. This resulted in yearly inventory-like maps (each Landsat pixel was assigned tree measurements from the nearest inventory plot), from which forest AGB and other attributes were calculated. Matasci, Hermosilla [51] developed a RF-based imputation model trained by six Lidar measurements of return height and percentage. Using this model, forest structural variables, including AGB, were then indirectly imputed across space and time.
Ma, Domke [80] introduced a novel approach using matrix models for estimating AGB dynamics from LTS. The study first developed a matrix model based on remeasured inventory data to predict AGB from forest structural variables including stand basal area. The basal area in the matrix was then replaced by Landsat spectral variables to project large-scale AGB dynamics short (5 years) and long-term (30 years). This modelling method incorporates various combinations of remote sensing and inventory data, avoiding a loss in accuracy from only using inventory variables. However, this approach is dependent on the availability of cyclical inventory data.

Accuracy Assessment
Depending on the availability of validation data, different approaches can be used to assess the performance of AGB models (Table 3). Where the sample size is large enough, studies often split the reference data into two sets (training and testing data), so the model can be cross-validated using the testing data. However, this approach often increases the variance in model training data, causing negative impacts on model accuracy. Therefore, where the sample size of reference data is small, studies often use internal assessment methods such as RF out-of-bag (OOB) errors and n-fold cross-validation (also known as leave-one-out cross-validation when n = 1). Some studies compare AGB perdition maps with independent AGB products such as national biomass maps [58,82] or Lidar-based AGB data [9,49,52]. Lidar-based AGB data are often used as they are widely confirmed as the most accurate remote sensing-based predictions.
Perhaps more important than model accuracy, when extending single model through time, is the need to evaluate (1) the temporal transferability of the model and (2) the robustness of the model in capturing AGB changes resulting from forest disturbance and recovery. Unfortunately, this is often dependent on the availability of validation data in multi-points in time (i.e., cycle inventory or Lidar data), which are normally not available in many forested regions. Few studies that estimated historical AGB have examined the temporal transferability of the model extension (Table 3). Powell, Cohen [81] performed a Landsat scene-level assessment by comparing scene-level predicted AGB trajectories with temporal trends of observed AGB derived from remeasured inventory plots. Matasci, Hermosilla [51] conducted a multi-temporal assessment of annual AGB predictions using Lidar data acquired from 2006 to 2012. In a recent study, Nguyen, Jones [52] conducted a time-series assessment of 30-year annual AGB prediction maps using un-changed, forested Lidar-based plots. Interestingly, these studies all confirmed that a single-date LTS-based model is transferable to estimate AGB over time.
Validating of AGB change is the most challenging as it requires repeated AGB observations. Where remeasured inventory plots are available, studies often extracted and compared observed AGB change (i.e., the difference of AGB between measurement rounds) with predicted AGB change [12,50,63,80,81,84]. Powell, Cohen [81] indicated that the accuracy of AGB change prediction improved when applying a pixel-based temporal fitting to AGB trajectories. Boisvenue, Smiley [50] fit a linear mixed-effect model to both observed and predicted AGB change datasets to estimate the long-term mean values of AGB change for each forest type (∆AGB = f (age)), which were then compared together. Where remeasured inventory plots are unavailable, multi-temporal Lidar data can be used as an alternative. Nguyen, Jones [52] did this when evaluating the ability of a LTS-based model in capturing AGB change resulting from forest disturbance and recovery processes. If repeated AGB observations are not available for validation, comparing with results from other studies could be an alternative. For example, Powell, Cohen [82] compared LTS-based AGB loss from disturbance (at a stratum-level) against results from several studies conducted within the same study area.

Characterising Spatial and Temporal Patterns of Forest AGB Dynamics
Forest AGB maps are normally created at the spatial resolution of 30 m, aligning with Landsat pixels, and at an annual or near-annual temporal density (few studies produced seasonal or multi-year epochal maps, Table 3). AGB change analysis across space and time is a crucial component of a forest biomass monitoring system. Along with developing robust modelling methods for estimating historical AGB, studies have proposed various approaches for quantifying spatial and temporal patterns of AGB dynamics to effectively support forest management (Table 3).
Methods for characterising patterns of AGB dynamics can be stratified into two broad groups: difference-based and trajectory-based. Difference is a traditional method for calculating AGB change by subtracting AGB values between two prediction dates. Though this is a simple and straightforward method, it often provides only generic information on AGB change. This method has been normally used when the temporal density of AGB predictions is low (often multiple-year epochal) and the intervals between prediction dates are wide enough to capture actual AGB changes in forests (5-10 years) [58,63,83,84]. As an exception, Boisvenue, Smiley [50] used the difference-based method to compute AGB changes from 29 annual AGB maps, but the authors then developed a mixed-effect model to summarize long-term trends of AGB at a stratum-level.
In contrast, most studies predicting AGB at an annual scale often used a trajectory-based method for characterizing AGB dynamics (Table 3). This method allows consistently and comprehensively tracking and understanding AGB changes across space and time according to forest disturbance and recovery processes. Some studies applied a temporal fitting process to annual AGB trajectories to create forest AGB disturbance and recovery maps [72,81]. Others integrated fitted AGB trajectories with LTS-derived forest disturbance maps to quantify AGB loss from disturbance [82] and according to different disturbance causal agents [9]. Matasci, Hermosilla [51] combined annual AGB predictions with an LTS-derived disturbance map to characterise AGB dynamics of Canadian forests across three decades. Temporal AGB trends for each ecozone were summarised for three scenarios: undisturbed forests, forests impacted by wildfires and harvesting. Nguyen, Jones [52] characterised spatial and temporal patterns of AGB dynamics by computing change metrics, including AGB loss and gain from disturbance and recovery, Recovery Indicator (RI), and Years to Recovery (Y2R). These allowed authors to robustly track AGB change according to different disturbance scenarios (i.e., disturbance severity and causal agents).

Conclusions and Future Opportunities
The opening to the public of the Landsat archive since 2008 has facilitated the use of LTS for forest applications. In this review, we provide an overview of recent studies utilizing LTS for estimating forest AGB and its dynamics. In summary, the use of LTS has the potential to improve the estimates of forest AGB in comparison with using single-date images. It provides a unique opportunity to create spatially and temporally complete products of forest AGB across large areas and long time-series. Furthermore, time-series processing and variables extracted from LTS can help improve the consistency and accuracy of AGB modelling. Many innovative approaches for estimating forest AGB across space and time have been demonstrated throughout the last decade, with methods becoming more advanced and robust in more recent studies. Thus, LTS-based products of AGB dynamics have become more reliable to support forest managers and researchers.
Estimating forest AGB dynamics across space and time is an ongoing topic of interest. LTS is increasingly used for this purpose as it provides the longest (and free) collection of imagery. Most current studies focused on boreal and temperate forests, highlighting the need for exploring the utility of LTS for estimating AGB dynamics in tropical forests. In contrast to boreal and temperate forests, using LTS to calculate AGB dynamics is all the more challenging in tropical forest regions due to persistent cloud cover. In addition, tropical forests are mainly found in developing countries, which often lack a systematic inventory network and Lidar data. Implementing a single model to predict AGB over time can produce temporal uncertainties (e.g., the balance of the model over time), which has mostly been ignored in current studies. More research investigating both spatial and temporal accuracies of LTS-derived AGB change are anticipated over coming years. Along with Landsat, other satellites such as Sentinel 2 and SPOT have been providing imagery as a time-series, with higher spatial and temporal resolutions. Incorporating such remote sensing data with LTS for improving the estimates of forest AGB dynamics will be a main area of research in the near future. Some studies are also investigating opportunities to combine LTS with active remote sensing data such as synthetic aperture radar (SAR) and the Global Ecosystem Dynamics Investigation (GEDI) for forest biomass estimation [99]. This is especially promising for tropical areas that suffer from extensive cloud/shadow issues. Last but not least, some scientists are attempting to project future forest AGB dynamics using LTS [80]. This is a challenging task for the forest remote sensing community, expected to be one of the main foci in the next few years. Trends of AGB over the last 40 years derived from LTS will be fundamental for predicting future AGB dynamics.

Conflicts of Interest:
The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; or the writing of the manuscript. Table A1. Landsat spectral indices commonly used for forest AGB estimates.

Landsat Spectral Index Calculation
Normalized Difference Vegetation Index (NDVI) [