Sub-annual tropical forest disturbance monitoring using harmonized Landsat and Sentinel-2 data

Accurate sub-annual detection of forest disturbance provides timely baseline information for understanding forest change and dynamics to support the development of sustainable forest management strategies. Traditionally, Landsat imagery was widely used to monitor forest disturbance, but the low temporal density of Landsat observations limits the timely detection of forest disturbance. Recently a new harmonized dataset of Landsat and Sentinel-2 imagery (HLS) has been created to increase the density of observations and provide more frequent images, but the added-value of this dataset for sub-annual tropical forest disturbance monitoring has not been investigated yet. Here, we used all available HLS images acquired from 2016 to 2019 to detect forest disturbance at two tropical forest sites in Tanzania and Brazil. Based on HLS data, forest disturbance was detected by combining normalized difference moisture index (NDMI) and normalized difference vegetation index (NDVI) time series using BFAST monitor and random forest algorithms. To assess the added-value of the HLS time series, we also detected forest disturbance from (i) Landsat-8/OLI time series only and (ii) Sentinel-2 time series only data. The spatial accuracy assessment of forest disturbance detection at the Tanzania site shows that the combined Landsat-8/OLI and Sentinel-2 data achieved the highest overall accuracy (84.5%), more than 3.5% higher than the accuracy of using only Landsat-8/OLI or Sentinel-2. Similarly, for the Brazil site, the overall accuracy of using the combined Landsat-8/OLI and Sentinel-2 data was 95.5%, at least 2% higher than others. In terms of temporal accuracy, the mean time lag of 2.0 months, was achieved from the combined data and Sentinel-2 only at the Tanzania site. This mean time lag is at least one month shorter than that of using Landsat-8/OLI only (3.3 months). At the Brazil site, the mean time lag of forest disturbance detection based on the combined data was 0.22 months, shorter by 0.50 and 0.15 months when compared to using Landsat-8/OLI only (0.72 months) or Sentinel-2 only (0.37 months), respectively. Our results indicate that HLS data is promising for accurate and timely forest disturbance detection particularly in the moist forest where cloud cover is often high.


Introduction
Characterizing forest disturbance at a sub-annual scale is of great importance for a better understanding of forest dynamics and assist in developing strategies for sustainable forest management. Forests cover approximately 30% of the land surface area and provide a broad range of ecological, economic and social services to human beings (Bonan, 2008). Forest disturbance disrupts the structure, composition and function of the forest ecosystem (Seidl et al., 2017) and it has been considered as an essential ecological process that is not well understood or quantified (Masek et al., 2008;Zhao et al., 2018). Detecting forest disturbance events at a sub-annual scale provides timely information for developing early intervention strategies to reduce illegal logging activities (Hamunyela et al., 2016a;Tang et al., 2019). For example, the National Institute for Space Research (INPE) of Brazil has developed an early warning system based on 56 m spatial resolution satellite image acquired from AWiFS (Advanced Wide Field Sensor) and 64 m spatial resolution imagery acquired from Wide Field Imager (WFI) sensor of China-Brazil Earth Resources Satellite 4 (CBERS-4), which is essential for preventing deforestation in the Brazilian Legal Amazon (Diniz et al., 2015). Similarly, a near-real-time forest disturbance detection also has been conducted in Congo Basin based on Sentinel-1 imagery (Reiche et al., 2021).
Satellite remote sensing data has been increasingly used in forest monitoring (Shimizu, Ota, and Mizoue, 2019). The capability of monitoring forest disturbances using satellite imagery with various spatial resolutions has been demonstrated in previous studies. For example, coarse-resolution imagery acquired from Advanced Very High Resolution Radiometer (AVHRR) or Moderate Resolution Imaging Spectroradiometer (MODIS) was adopted to detect forest disturbances (Jin and Sader, 2005;Hammer et al., 2014). Similarly, images obtained from Landsat and Sentinel-2 have also been successfully used to monitor forest disturbances and regrowth (Potapov et al., 2015;Wang et al., White et al., 2018;Lima et al., 2019). Also, very high-resolution imagery (VHR), like IKONOS has been used to characterize forest disturbances (Frolking et al., 2009). In addition, VHR imagery (e.g. WorldView-2 and GeoEye-1) has been combined with LiDAR to detect canopy tree loss in the Jamari National Forest, Brazil (Dalagnol et al., 2019).
Landsat images were frequently used in forest monitoring due to their free availability and long-term monitoring capability (Wulder et al., 2019). Many attempts have been made to monitor forest disturbances at annual or sub-annual scales using Landsat time-series data (Kennedy et al., 2010;White et al., 2017;Zhao et al., 2018). For example, Kennedy et al. (2012) have studied the disturbances and recovery in the United States based on Landsat time-series data. Similarly, forest disturbances and regrowth in Canada have been characterized annually from 1985 to 2010 using vegetation indices derived from the Landsat time series (White et al., 2017). In the Brazilian State of Rondônia, forest disturbances were recently detected using spectral unmixing and Landsat time-series data (Bullock et al., 2020). Although Landsat images are widely used for forest monitoring, the low temporal density of Landsat observations limits the sub-annual detection of tropical forest disturbance. Tropical forests are located in the cloud persistent regions, which makes the actual temporal revisit cycle of Landsat-8/OLI images sometimes longer than 16 days (Shao et al., 2019). With the launching of Sentinel-2A and Sentinel-2B sensors in June 2015 and March 2017, satellite images with 10 m to 60 m spatial resolution and a revisit time of 5 days at the equator are freely accessible (Drusch et al., 2012). Integrating imagery from Sentinel-2 sensors with Landsat imagery would make denser time series data available (Claverie et al., 2018;Shao et al., 2019).
Previous studies have suggested that combining Landsat and Sentinel-2 imagery has multiple advantages for land monitoring. Combining Landsat-8/OLI and Sentinel-2 imagery can provide more frequent data coverage. For example, when combining Landsat-8/OLI, Sentinel-2A and Sentinel-2B (three sensors), the maximum revisit interval is approximately 7.0 days (Li and Roy, 2017). Furthermore, combining Landsat and Sentinel-2 is promising for characterizing changes or dynamics with more details . Dense time series data has the capability of providing more detailed information compared with bi-temporal or annual satellite data and it is suitable for monitoring forest disturbance in the tropics as the vegetations recover rapidly after disturbance (Shimizu et al., 2019). Harmonized Landsat Sentinel-2 (HLS) data was developed by a National Aeronautics and Space Administration (NASA) project that aims to provide seamless product from surface reflectance data acquired by Landsat-8/OLI imagery (acquired since 2013) and Sentinel-2 imagery (acquired since October 2015) (Claverie et al., 2018). HLS data consists of three products, S10 (10 m, 20 m and 60 m spatial resolution) derived from Sentinel-2 imagery, S30 (30 m spatial resolution), derived from S10 product and L30 (30 m spatial resolution) derived from Landsat-8/OLI. While the HLS has been successfully used to monitor land surface phenology and characterize grassland use intensity (Bolton et al., 2020;Griffiths et al., 2020), its potential for forest disturbance monitoring has not been examined yet.
Numerous change detection algorithms are widely used in detecting forest disturbances. Landsat-based detection of Trends in Disturbance and Recovery (LandTrendr) is a temporal segmentation-based approach that extracts spectral trajectories of land surface from yearly Landsat time-series data (Kennedy et al., 2010). Vegetation Change Tracker (VCT) is an automated approach that takes advantage of the spectraltemporal properties of the Landsat time series stack to reconstruct the disturbance history of forests (Huang et al., 2010). Continuous Change Detection and Classification (CCDC) uses all available Landsat data to characterize the land cover change in near-real-time (Zhu and Woodcock, 2014). The Breaks For Additive Seasonal and Trend (BFAST) proposed by Verbesselt et al. (2010), decomposes time series into trend, seasonal and remainder to detect structural changes in any time series data. BFAST Monitor, developed by Verbesselt et al. (2012) is used to detect disturbances in near real-time monitoring while the LandTrendr and VCT produce annual maps. Both CCDC and BFAST Monitor can be used for sub-annual forest disturbance. BFAST Monitor was selected in this study because it has been successfully used to detect small-scale forest disturbances in the tropical montane forest (DeVries et al., 2015).
This study aims to investigate the added value of integrating Sentinel-2 and Landsat time-series data for sub-annual tropical forest disturbance monitoring. More specifically, we address the following questions: (1) What is the spatial and temporal accuracy of forest disturbance detected from Landsat-8/OLI, Sentinel-2 or the combined data? (2) Whether the use of the combined Landsat-8/OLI and Sentinel-2 data improves the spatial and temporal accuracy of forest disturbance detection compared with using only single sensor data or not? In this paper, forest disturbance is defined as the complete removal of tree cover within a 30 m by 30 m pixel (~0.1 ha). We followed the same definition that was used in other satellite-derived products of forest change monitor (Hansen et al 2013;Vargas et al., 2019;Reiche et al. 2021). We used BFAST Monitor for structural change detection and combined it with random forest to predict the forest disturbance in two study sites located in the tropical regions (Brazil and Tanzania). We tested the accuracies of forest disturbances detection using L30 (derived from Landsat-8/OLI), S30 (derived from Sentinel-2) and the combination of both products.

Study area
The study area consists of two HLS tiles (same as the tile of Sentinel-2 L1C data), 37MCP (covers an area of 11,938 km 2 ) and 21LYH (covers 12,018 km 2 ), located in Tanzania and Brazil, respectively ( Fig. 1). Tile 37MCP is located at the intersection of Morogoro, Pwani and Tanga regions in Tanzania. Tile 21LYH is located within the State of Mato Grosso, Brazil. According to the Terrestrial Ecoregions of the World, the Brazil site (Tile 21LYH) is covered by seasonal forest and moist forest (Olson et al., 2001), whereas the Tanzania site (Tile 37MCP) is mainly covered by tropical try forest. Tropical dry forest is defined as vegetation dominated by deciduous trees that grow in places where the mean monthly temperature is ≥ 25•C, annual precipitation ranges from 700 mm to 2,000 mm and with three or more dry months (Sánchez-Azofeifa et al., 2005). The HLS tile 37MCP has been selected because the montane forest in the Mvomero district (one district of the Morogoro region) has been disturbed by 27% from 2001 to 2017 (Hamunyela et al., 2020). Tile 21LYH in the State of Mato Grosso was selected because Mato Grosso is one of the deforestation frontiers in Amazon in recent decades (Gollnow et al., 2018). In the Mvomero district, the forest loss is driven by the expansion of commercial cropland, small-holder mixed crop-livestock farming and smallholder farming (Hamunyela et al., 2020). Disturbed forest in Mato Grosso is mainly converted into a large-scale pasture and agricultural land (Barona et al., 2010;Wang et al., 2019).

HLS data
NASA provides the HLS version 1.4 data, harmonized Landsat-8/OLI and Sentinel-2 data through a series of processing steps include radiometric adjustments, geometric adjustments, atmospheric correction, cloud, cloud-shadow masking and spatial co-registration, Bidirectional Reflectance Distribution Function (BRDF) normalization and bandpass adjustment (Claverie et al., 2018). The HLS data includes Landsat-8/OLI and Sentinel-2 imagery acquired since April 2013 and October 2015, respectively. HLS data consists of S10, S30 and L30 products. S10 product is not used in this study because its spatial resolution is different from other products. All available HLS images of the two HLS tiles (37MCP, 21LYH) from 2016 to 2019 were used in this study. Overall, 528 scenes (89 from Landsat-8/OLI and 439 from Sentinel-2) of images were used in the Tanzania site, and. 449 scenes (90 from Landsat-8/OLI and 359 from Sentinel-2) for the Brazil site.

Ancillary data
Ancillary data that were used in this study includes global forest change data (version 1.7) provided by Hansen et al. (2013) and Plan-etScope imagery (https://www.planet.com/). Global forest change data was used for delineating the benchmark forest mask. PlanetScope imagery was used for collecting training and test data for random forest models to predict the forest disturbance. PlanetScope imagery acquired by PlanetScope satellite constellation consists of more than 120 sensors (Wang et al., 2020) and has daily revisit time at nadir. PlanetScope imagery consists of four bands, including Blue (455-515 nm), Green (500-590 nm), Red (590-670 nm) and NIR (780-860 nm) bands (Cooley et al. 2017). All bands were used for interpreting images when collecting training and test data for the random forest model. The spatial resolution of PlanetScope imagery can achieve 3.7 m.

Methods
The workflow of detecting forest disturbance using the HLS is shown in Fig. 2. First, we pre-processed the HLS satellite images (see Section 2.4.1). The pre-process step mainly refers to the masking of cloud and cloud shadow. Secondly, a benchmark forest mask was generated to mask out the non-forest pixels. Next, two indices, Normalized Difference Moisture Index (NDMI) (Gao, 1996) and Normalized Difference Vegetation Index (NDVI) (Rouse et al., 1974) were calculated and stacked. Then, BFAST Monitor and random forest were applied to detect the forest disturbance (see Sections 2.4.2 and 2.4.3). Finally, the accuracies of forest disturbance maps based on different data were assessed (see Section 2.4.4).

Image pre-processing
2.4.1.1. Screening cloud and cloud shadow. HLS contains a quality assurance (QA) layer that identified the cloud, cloud shadow, snow and water generated from the Fmask algorithm. QA layer in each image was used to mask the cloud, cloud shadow, water and snow for the L30 product. However, the Sentinel-2 QA layer contains errors in cloud and cloud shadow detection (Bolton et al., 2020). For the S30 product, Fmask 4.0 developed by (Qiu et al., 2018) was applied to mask the cloud and cloud shadow. In addition, a threshold-based approach was used to further mask out the remaining cloud whereby a difference between NDVI and NDMI, larger than − 0.4 was considered as the cloud. The remaining shadows were masked using the shadow index (Zhu and Helmer 2018) whereby values were smaller than 0.34 were considered as cloud shadow. A trial-and-error approach was used to select thresholds for masking cloud and cloud shadows. Images with high cloud coverage were excluded for further analysis.

Benchmark forest mask.
An initial forest mask was produced from a random forest classification based on the composite imagery of 2018. The composite imagery is based on the median value of all available Sentinel-2 imagery acquired from June to December 2018. Both study sites are located in cloud persistent tropical regions and almost half a year's data were needed to cover the whole study area. Sentinel-2 images were used, instead of Landsat-8/OLI images, to generate the forest mask because they have red edge bands that can provide additional information on forest cover. A total of 420 training sampling pixels (210 for the forest, 210 for non-forest) were used to train a random forest model for mapping the forest. Furthermore, the Global Forest Change data was used to generate the second forest mask. Tree canopy cover for the year 2000 and the forest cover loss and gain data during the period 2000-2018 were also utilized. The tree cover that is less than 25%, forest loss during 2000-2018 was masked out, forest gain during 2000-2018 was added. Finally, to reduce the committed error of forest class in the map produced from the image classification, we further refined the forest mask by masking out the pixels that were identified as forest in the initial forest mask whereas labeled as nonforest in the second forest mask.

Structural change detection based on BFAST Monitor
BFAST Monitor was used to detect the structural changes in the time series of NDMI and NDVI for Landsat-8/OLI only, Sentinel-2 only and the combined data. BFAST Monitor was selected in this study because it has been successfully used to detect small-scale forest disturbances in the tropical montane forest (DeVries et al., 2015). To apply the BFAST Monitor, we used "bfastSpatial", a wrapper for the BFAST monitor. The implementation of BFAST Monitor consists of several steps, season-trend model fitting, structural change detection, the calculation of Magnitude and R 2 .
Each pixel time series was divided into history period and monitoring period (Fig. 3). The history period in this study was between the year 2016 and 2018 and the monitoring period was the year 2019. We then applied a season-trend model (BFAST Monitor) described in Verbesselt et al. (2012) in the history period. For the observation y t at time t, a season-trend model is fitted: where α 1 is the intercept, α 2 is the slope (trend), γ j and δ j are amplitudes and phases,k is the chosen harmonic order, f is the temporal frequency and ε t is the random error (Verbesselt et al., 2012). We used the firstorder model and omitted the trend component in the formula because including the trend may lead to false break detection and inflate magnitude (DeVries et al., 2015). The model can be written as a standard linear regression model and more details were explained in Verbesselt et al. (2012). Ordinary Least Squares (OLS) was used to estimate coefficients in the linear model based on the observations (Zeileis et al., 2002).
Based on the described model, assuming that the model is stable for an observed time period, if the new observations make the model unstable, a structure change is detected. For structure change detection, in the monitoring period, moving sums of residuals were calculated by: where σ is the estimator of the variance, n represents the number of observations in history period, h denotes the bandwidth that determines the moving data window size, y s and ŷ s represent the actual observations and expected observations, respectively (DeVries et al., 2015;Verbesselt et al., 2012). A breakpoint was detected if the absolute value of MO t is significantly different from zero at a specified confidence level (DeVries et al., 2015;Verbesselt et al., 2012;Zeileis et al., 2002). Magnitude value is the median of the difference between the observed value and predicted value for all the observations in the monitoring period. For each pixel (with detected breaks or without detected breaks), the magnitude value was calculated. The median value instead of the sum value was calculated because the median is less sensitive to noise (DeVries et al., 2015). R 2 is derived from the linear regression model used in the history period.

Random forest for detecting forest disturbance
BFAST Monitor can return the breakpoint, magnitude and R 2 . The magnitude is the median of the difference between the observed value and predicted value for all the observations in the monitoring period. Using the BFAST Monitor approach for forest disturbance detection, an additional filtering approach needs to be conducted to extract forest disturbance (Hislop et al., 2019). The breakpoint and magnitude are usually post-processed using either a threshold-based approach or random forest for identifying possible forest disturbances (DeVries et al., 2015;Schultz et al., 2018). Previous research has tested these two approaches and found out that the random forest algorithm outperformed the threshold-based approach (Schultz et al., 2018). Therefore, a random forest algorithm (Breiman, 2001) was used to predict the forest disturbance based on the break detection results from NDMI and NDVI time series. Random forest classification was used to differentiate the forest disturbance from the undisturbed forest. For the random forest model, the input features include six raster layers, the presence/absence of breakpoint (a binary variable converted from the breakpoint layer), magnitude value and R 2 from the BFAST Monitor output based on NDMI and NDVI time series. A total of 1000 pixels (500 for each site) were selected and interpreted based on the original images and the Planet-Scope imagery. Among the 500 pixels for each site, 250 were forest disturbance pixels and the remaining were undisturbed pixels. For the random forest model, the number of decision trees was set to 500. The number of randomly selected variables when splitting at the tree node was set to 2 (the square root of the number of input features). Random forest models were implemented using the R package "randomForest" (Liaw and Wiener, 2002) for three data products of both study sites.

Accuracy assessment
The spatial and temporal accuracies of the forest disturbance detected in 2019 were assessed. For the spatial accuracy, the disturbance maps were validated based on stratified random sampled pixels (30 m by 30 m) created based on the forest disturbance maps, generated from classification using Landsat-8/OLI, Sentinel-2 and the combined data. 400 random pixels were selected for each disturbance map. Each forest disturbance map consists of two classes (forest and forest disturbance), for each class, 200 random pixels were generated and interpreted. The occurrence of forest disturbance in 2019 was interpreted based on PlanetScope imagery. 1200 sample pixels for each site and 2400 sample pixels for the total of two sites. Based on the confusion matrix, commonly used producer's accuracy, user's accuracy and overall accuracy, described by Congalton (1991), were calculated. For assessing the temporal accuracy, the same validation pixels used to assess the spatial accuracy were used. The sampling points that were identified as forest disturbance in the map and in the validation data were selected for further analysis. The monthly mosaic of PlanetScope imagery was used to label which month the forest disturbance happened. When the BFAST monitor returned two dates of breakpoint detected based on NDMI or NDVI, the dates from NDMI and NDVI were compared and the earlier date was considered as the disturbance time. When the break was only detected from one index, NDMI or NDVI, the returned date was considered as the disturbance time. Then, the daily date was aggregated into monthly data in R. The disturbance time detected by the model and the validation data were compared. If the disturbance detected in the model was later than the validation data, the sampling pixel was considered as correctly identified. We also considered those disturbances detected in the model that are one month earlier than the reference disturbance date as correct detection. The mean time lag is defined as the mean days between the reference date and the date predicted by the model (Shimizu et al., 2019). In this study, we calculated the mean time lag between the model predicted time and reference data based on the month when the disturbance occurred or was detected by the model.

Forest disturbance detection
Forest disturbance detected from two study sites is shown in Fig. 4. In the Tanzania site, the detected forest disturbance area using Sentinel-2 imagery was the largest (4,387 ha), followed by using the combined data (3,712 ha) and using only Landsat-8/OLI imagery (2,955 ha). Most of the forest disturbance was detected in the eastern part of the Tanzania site using Sentinel-2 and the combined data, while the detected forest disturbance was distributed in the eastern and southwest part using Landsat-8/OLI imagery. While for the Brazil site, the area of the detected disturbance using Landsat-8/OLI, Sentinel-2 and the combined data were 9,770 ha, 10,090 ha, and 8,706 ha, respectively. The spatial distribution of the detected forest disturbance using Landsat-8/OLI, Sentinel-2 and the combined data were similar. A large amount of forest disturbance was detected in the south part of the study area. Most of the detected forest disturbance sizes detected in the Brazil site were

Accuracy of the forest disturbance maps
The overall spatial accuracy of forest disturbances in the Tanzania site using the combined Landsat-8/OLI and Sentinel-2 is the highest (84.5%), followed by using only Landsat-8/OLI (81.0%) and Sentinel-2 (79.8%). In terms of the producer's accuracy of the disturbance class in the Tanzania site, the highest accuracy (87.1%) was achieved by using the combined data. Thus the omission error of the disturbance class is the smallest when using the combined data. For the forest disturbance class, using only Landsat-8/OLI data achieved the highest user's accuracy (84.0%), indicating the smallest commission error of disturbance class in the Tanzania site. The user's and producer's accuracies of the disturbance map created by using the combined data were more balanced compared to maps created based on Landsat-8/OLI or Sentinel-2 imagery in the Tanzania site.
The overall spatial accuracy of forest disturbances in the Brazil site using the combined Landsat-8/OLI and Sentinel-2 was the highest (95.5%), followed by using only Landsat-8/OLI (93.3%) and Sentinel-2 (92.0%) ( Table 1). The forest disturbance map based on the combined data achieved the highest overall accuracy, followed by using Landsat-8/OLI and using Sentinel-2. In terms of the producer's accuracy of the disturbance maps, the highest accuracy of change class (disturbance) was 98.8%, achieved by using only Sentinel-2. The highest user's accuracy of the disturbance class was 92.5%, which was achieved by using the combined data. For the disturbance class, the commission errors were relatively higher than the omission errors.
Among the sampled pixels for validating the temporal accuracy in the Tanzania site, the number of pixels that was correctly identified by using Landsat-8/OLI and Sentinel-2 were 151 (89.9%) and 117 (82.4%), respectively. While for the combined data, the number was 140 (86.4%). The mean time lags for forest disturbance mapping using Landsat-8/OLI and Sentinel-2 was approximately 2.0 months (Fig. 6). When detecting the forest disturbance using Landsat-8/OLI, the mean time lag was 3.3 months. The mean time lag for forest disturbance detection using Landsat-8/OLI was the longest.
Regarding the temporal accuracies in the Brazil site, the number of pixels that was correctly identified by using Landsat-8/OLI, Sentinel-2 and the combined data were 164, 136, 151, respectively and they accounted for 93.2%, 80.0% and 81.6% of the total sampled pixels used in the temporal accuracy assessment. The mean time lags for forest disturbance mapping using Landsat-8/OLI, Sentinel-2 and the combined data were approximately 0.72 months, 0.37 months, and 0.22 months, respectively (Fig. 6). When using the combined data for forest disturbance detection, the mean time lag was smaller compared with disturbances detected by using Landsat-8/OLI or Sentinel-2 imagery.

Discussion
The added-value of using HLS data to detect disturbances in tropical forest systems at a sub-annual scale was assessed based on BFAST Monitor and random forest models for two study sites located in Tanzania and Brazil. Our results show that using the combined data achieves a shorter mean time lag when detecting forest disturbance compared with using only Landsat-8/OLI. Furthermore, combining the Landsat-8/OLI and Sentinel-2 imagery also improves the overall spatial accuracy for detecting forest disturbance.

Spatial accuracy
Detecting forest disturbance using the combined Landsat-8/OLI and Sentinel-2 imagery produced the highest spatial accuracy. In both study sites, the spatial accuracy of forest disturbance detection using the  combined data was higher than using only Landsat-8/OLI or Sentinel-2 data. The higher overall accuracy of forest disturbance detection using the combined data can be attributed to the increased number of valid observations in the pixel time series by integrating imagery from different sensors (Fig. 7A). A larger number of valid observations lead to a better fit of the BFAST Monitor model thus achieving higher accuracies (Schultz et al., 2016). Previous research compared the performance for using Sentinel-2 and Landsat-8/OLI imagery to monitor selective  logging in the Brazilian Amazon and found out that the overall accuracy of the disturbance and undisturbed forest map based on Sentinel-2 was marginally higher than the one based on Landsat-8/OLI (Lima et al., 2019). On the contrary, our results show that forest disturbance detection using Landsat-8/OLI achieves higher accuracy compared with using Sentinel-2 data. This can be attributed to two reasons. On the one hand, the original Sentinel-2 imagery was resampled to a 30 m spatial resolution imagery in the HLS product, so the advantage of Sentinel-2 in revealing better spatial details was lost. On the other hand, unmasked cloud-contaminated pixels in Sentinel-2 imagery can affect forest disturbance detection. Fmask 4.0 was adopted to mask out the cloud, but Fmask 4.0 produces better results in detecting cloud by using Landsat-8/ OLI compared with using Sentinel-2 because Sentinel-2 does not include the thermal bands (Qiu et al., 2019).

Temporal accuracy
Regarding the temporal accuracy of the forest disturbance detection, combining Landsat-8/OLI and Sentinel-2 data produced a shorter mean time lag than using Landsat-8/OLI imagery. The mean time lag decreased when using the combined data because of the increased temporal resolution compared with using only Landsat-8/OLI imagery. Similarly, a forest disturbance detection study using the combined Sentinel-1, ALOS-2 PALSAR-2 and Landsat data showed that the temporal accuracies based on multi-sensors were higher than single-sensor (Reiche et al., 2018). Our findings show that integrating Landsat-8/ OLI and Sentinel-2 enables earlier forest disturbance detection. The value of using dense time-series data (HLS) for forest disturbance monitoring in the tropics is apparent and it can contribute to building an early warning system.
The mean time lags of forest disturbance in the Brazil site are relatively shorter compared with those in the Tanzania site, which might be explained by the difference in the number of valid observations, patch sizes of forest disturbances and forest types. For each sensor, the number of valid observations of Landsat-8/OLI, Sentinel-2 and the combined data in the Brazil site is larger than those in the Tanzania site (Fig. 7A). Observation availability is the driving factor for evaluating the performance of forest disturbance monitoring systems based on satellite images (Reiche et al., 2018). Future studies can investigate the potential of adding the Landsat time-series data acquired from 2013 to 2015 to increase the number of valid observations in sub-annual forest disturbance monitoring. In addition, the patch sizes of forest disturbances might affect the time difference of forest disturbance detection in different sites. The number of large patches of forest disturbances in the Brazil site is higher than that of the Tanzania site (Fig. 7B). The seasonal variation of forests in the Tanzania site is larger than that in the Brazil site (Fig. 3). Future research can investigate the use of a spatial normalizer to reduce the seasonal variations (Hamunyela et al., 2016a(Hamunyela et al., , 2016b. The temporal accuracy assessed in this study was limited by the validation data. When assessing the temporal accuracy of detecting forest disturbance using monthly mosaic PlanetScope imagery as the reference data may introduce uncertainty because there is a discrepancy between the time when the disturbance happens and the time when a satellite image captures the disturbance. To reduce the uncertainty, acquiring satellite images at a higher temporal resolution is necessary for validating the sub-annual detection of forest disturbance. With the launching of more satellites, more frequent observations are available, which will increase the observation availability of the reference images (Shimizu et al., 2019). Besides, an adjusted reference approach can be applied to reduce the imprecision caused by the time lag between the change truly occurs and the change captured by the images (Reiche et al., 2018).
Visual comparison with the PRODES product in the Brazil site shows that the forest disturbances based on the combined Landsat-8/OLI and Sentinel-2 data provides complementary information for the PRODES. PRODES, one of the main monitoring projects conducted by INPE, updates annual deforestation maps in the Legal Amazon since 1988 using largely Landsat imagery (Hansen et al., 2008). There is an overall agreement in the distribution of detected disturbed forests ( Fig. 8A and B). The major differences (blue rectangles in Fig. 8) are the disturbances detected only from PRODES and the difference could be related to the use of different benchmark forest masks. PRODES updates its annual map from August 1, while our forest disturbance detection started from January 1. Most of the detected disturbances within blue rectangles of PRODES ( Fig. 8C and D) could have occurred before 2019. Furthermore, forest disturbances detection using HLS data and BFAST Monitor approach has the capability of providing monthly forest disturbance products, which reveals more detailed information on when the forest disturbance occurs compared with the annual disturbance product of the PRODES. Detecting forest disturbance using the combined Landsat-8/ OLI and Sentinel-2 data not only has the potential of providing complementary information on which month the forest disturbance occurred but also enables earlier detection of forest disturbance (Fig. 6). Detecting the forest disturbance as early as possible provides significant information for policymakers to rapidly respond to reduce and prevent the expansion of forest disturbance.

Major sources of errors
Clouds, cloud shadows and forest masks affect the performance of the forest disturbance detection approach used in this study. For unmasked cloud and cloud shadows, a single outlier can be handled by BFAST Monitor and it is unlikely to be identified as a breakpoint although it decreases the NDMI/NDVI value because the approach using the moving sums of residuals and the calculation of magnitude value is based on the median value within the monitoring period (DeVries et al., 2015). However, unmasked cloud-contaminated or shadowcontaminated pixels that occur in multiple continuous images challenge the breakpoint detection. Integrating imagery acquired from other optical or Synthetic Aperture Radar (SAR) sensors may reduce the number of persistent cloud or shadow contaminated pixels. An accurate benchmark forest mask plays an important role in forest disturbance monitoring (Hirschmugl et al., 2020). If the benchmark forest mask does not exclude any forest disturbances that occurred before the observation period (e.g., 2019), these pixels would be mislabeled as a disturbance at the beginning of the monitoring period, a similar issue has been found in previous research (DeVries et al., 2015). Analyzing the stability of the data in the history period and using the stable observations instead of all the observations to fit the season-trend (Verbesselt et al., 2012) might reduce the error caused by inaccurate forest mask.

Implications for forest disturbance monitoring
Using HLS data to detect forest disturbance at a sub-annual scale would make a significant contribution to the disturbance monitoring system, particularly in terms of identifying the time of forest disturbance. For example, monitoring projects such as the PRODES project could benefit from the combined satellite data. The PRODES project updates its deforestation map annually. While our results produce subannual products for forest disturbance monitoring, which provide complementary information on which month the forest disturbance occurs. The average mean time lags of forest disturbance detection in the Brazil site are less than 1 month. This indicates that our results have the potential of providing information for rapid response on preventing further carbon or biodiversity loss. Furthermore, the minimum map unit of PRODES is 6.25 ha, a large proportion of the number of forest disturbances events that are less than 1 ha were detected in the Brazil site in 2019 (Fig. 7B) using the combined Landsat and Sentinel-2 imagery, which also provides complementary information for PRODES. Other Landsat-based disturbance alert system, such as the humid tropical forest disturbance monitoring system implemented in Peru, the Republic of Congo and Kalimantan, Indonesia (Hansen et al., 2016) would benefit from the combined data because integrating Sentinel-2 imagery and Landsat-8 OLI imagery for forest disturbance monitoring enables earlier forest disturbance detection. The earlier the disturbance is detected, the earlier the intervention action could be taken.
For using HLS imagery to monitor forest disturbance in other regions, several key issues needed to be addressed. First and foremost, the added value of using the combined Landsat-8/OLI and Sentinel-2 data for forest disturbance detection is more pronounced in the temporal component compared with the spatial accuracy in the two test study sites. Secondly, although the combined Landsat-8 OLI and Sentinel-2 data increase the temporal resolution, it also decreases the spatial resolution of Sentinel-2 imagery, which makes it inappropriate for detecting small-scale size forest disturbance that is less than one pixel (30 m by 30 m). For small-scale forest disturbance, 10 m resolution Sentinel-2 imagery or other higher spatial resolution satellite data are recommended. Future research could investigate the trade-offs introduced by HLS between the gain in temporal resolution and loss of spatial resolution (Sentinel-2 imagery) in forest disturbance monitoring, which is particularly important for small-scale forest disturbance detection in smallholder-dominated landscapes in Africa, such as Tanzania. Thirdly, for upscaling HLS imagery to monitor forest disturbance at a larger scale (e. g. national or global scale), a large amount of representative training and test samples are needed to be selected and interpreted. Different data availability and heterogeneity of landscape also pose challenges for a larger scale forest disturbance monitoring. Lastly, the use of HLS data for forest disturbance monitoring has its advantages and disadvantages compared with SAR-based or combined SAR and optical imagery based disturbance monitoring. The main disadvantage of using HLS imagery is the cloud contamination problem. Regarding the advantages, the HLS does not have speckle in the image and it usually has higher temporal resolution compared with most SAR imagery. Although integrating HLS with SAR improves the temporal resolution, the inconsistent data availability, distinct physical properties between SAR and optical images might be challenging (Markert et al., 2018).

Conclusions
Understanding forest disturbance provides baseline information for quantifying the carbon balance of the forest ecosystem and is critical for developing sustainable forest management strategies. We investigated the potential of using dense time-series HLS data for sub-annual tropical forest disturbance detection in two different sites. Our results show that integrating Sentinel-2 imagery with Landsat-8/OLI imagery has the potential of improving forest disturbance monitoring. Combining the Landsat-8/OLI and Sentinel-2 imagery can detect the forest disturbance with marginally higher spatial accuracy compared with using only Landsat-8/OLI imagery or using only Sentinel-2 imagery. More importantly, combining the Landsat-8/OLI and Sentinel-2 imagery shortens the mean time lag in forest disturbance monitoring. With more dense time-series data acquired from Landsat 10, Sentinel-1, PlanetScope or other sensors, the accuracy of sub-annual forest disturbance detection is expected to improve further and computing power needs to be increased to combine data from multiple sensors. HLS data is therefore promising for tropical forest disturbance monitoring and it opens up a new opportunity for forest disturbance detection in other study sites when global coverage of HLS Version 1.5 data are available.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. constructive comments on the manuscript. We also acknowledge NASA for the provision of the HLS data.