Mapping Forest Stability within Major Biomes Using Canopy Indices Derived from MODIS Time Series

: Deforestation and forest degradation from human land use, including primary forest loss, are of growing concern. The conservation of old-growth and other forests with important environmental values is central to many international initiatives aimed at protecting biodiversity, mitigating climate change impacts, and supporting sustainable livelihoods. Current remote-sensing products largely focus on deforestation rather than forest degradation and are dependent on machine learning, calibrated with extensive ﬁeld measurements. To help address this, we developed a novel approach for mapping forest ecosystem stability, deﬁned in terms of constancy, which is a key characteristic of long-undisturbed (including primary) forests. Our approach categorizes forests into stability classes based on satellite-data time series related to plant water–carbon relationships. Speciﬁcally, we used long-term dynamics of the fraction of photosynthetically active radiation intercepted by the canopy (fPAR) and shortwave infrared water stress index (SIWSI) derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) for the period 2003–2018. We calculated a set of variables from annual time series of fPAR and SIWSI for representative forest regions at opposite ends of Earth’s climatic and latitudinal gradients: boreal forests of Siberia (southern taiga, Russia) and tropical rainforests of the Amazon basin (Kayap ó territory, Brazil). Independent validation drew upon high-resolution Landsat imagery and forest cover change data. The results indicate that the proposed approach is accurate and applicable across forest biomes and, thereby, provides a timely and transferrable method to aid in the identiﬁcation and conservation of stable forests. Information on the location of less stable forests is equally relevant for ecological restoration, reforestation, and proforestation activities.


Introduction
Deforestation and forest degradation from anthropogenic impacts are now understood as a matter of global environmental concern [1]. Primary forests, which represent the most natural end of the forest ecological condition gradient, are under particular threat. According to the Food and Agriculture Organization (FAO) and the International Union for Conservation of Nature (IUCN), primary forests are defined as naturally regenerated forests of native species where there are no clearly visible indications of modern human land use activities and ecological processes dominate [2,3]. Primary forests typically have larger ecosystem carbon stocks, provide habitat for a wider range of species, yield cleaner water, and are overall more stable and resilient to external stressor, as well as have higher natural adaptive capacity, compared to degraded and secondary growth forest [4][5][6][7][8][9][10][11]. Therefore, the degradation of primary forests due to human land use and related impacts can result in Here we examine the constancy property of stability in tropical and boreal forests by using a novel approach based on remote-sensing indices that are known to be correlated with ecophysiological processes driving plant water-carbon relationships [33], which, in turn, are related to vegetation structure [34]. Specifically, we use the fraction of photosynthetically active radiation intercepted by the sunlit canopy (fPAR) [35,36] and the shortwave infrared water stress index (SIWSI) derived from the Moderate Resolution Imaging Spectroradiometer (MODIS). Moreover, fPAR is related to both vegetation structure (leaf area and orientation) and primary productivity [37,38], and SIWSI is a robust indicator of plant water content in fine tissues [39], particularly in green leaves [40]. SIWSI time series can therefore be indicative of water stress dynamics and the microclimate buffering capacity of structurally intact forest canopies. The latter determines the extent to which below-canopy climatic conditions deviate from the climate outside forests. We developed a set of variables calculated from annual time series of fPAR and SIWSI for two representative forested regions at opposite ends of Earth's climatic and latitudinal gradient: the boreal forests of Siberia (southern taiga, Russia) and tropical rainforests of the Amazon basin (Kayapó territory, Brazil). We validated our results by comparing them with several independent data sources (see Section 2.3 for more information) to ensure that our mapped stability categories equated to known distributions of disturbance and intact forest canopies (noting the inherent limitation of validating low stability forests, for which our method provides a first estimate).

Study Regions
We focused on two case study regions in contrasting forest biomes from different continents: the boreal taiga (Eurasia) and the tropical rainforest (South America). Within these, we targeted the southern taiga ecoregion of Siberia (Russia) and the Kayapó indigenous territory within the Amazon basin (Brazil) (Figure 1). These regions are representative of the aforementioned biomes and contain the full spectrum of forest ecosystem conditions: significant areas of primary forests [3] and intact forest landscapes still exist [41], along with forests that are heavily managed for commodity production and forests cleared for agricultural use. Many areas within both regions have been heavily exploited since the mid-20th century, resulting in substantial fragmentation and degradation. Fires and intensifying industrial logging are major threats to the sustainability of these forests, with the majority of fire ignitions in both regions caused by humans [42,43]. Moreover, human activities are radically altering their natural fire dynamics owing to altered fuel loads, desiccation, and flammability within or in proximity to human-altered landscapes [44][45][46][47][48][49][50].
Climate differs markedly between these regions, particularly with respect to thermal regime, degree of seasonality and growing season length (Supplementary Table S1). Extreme continentality is characteristic of the southern taiga, with large daily and seasonal ranges of temperature, solar radiation, and day length. These fluctuations are far smaller in tropical forests. Such differences have profound implications for forest composition, structure, function, and dynamics.
For the southern taiga, we defined the study area based on ecoregion boundaries to facilitate the identification of forest stability classes relative to a broad ecosystem type. In this way, we controlled for the existing differences in vegetation structure and floristics (e.g., dominant species, maximum canopy height and density) across boreal Russia, since these are largely influenced by the predominant physical environmental conditions. We used an ecoregion division for Russia from previous studies [51,52] that includes-from west to east-the administrative territories of Sverdlovsk, Tyumen, Omsk, Novosibirsk and Tomsk oblasts, Krasnoyarsk krai, and Irkutsk oblast. For the Amazon rainforest, the study region coincided with the Kayapó indigenous territory in the Brazilian states of Pará and Mato Grosso. This is an extensive area of largely primary tropical forest. We added a buffer zone of 60 km around the perimeter of the Kayapó territory to include areas that have been For the southern taiga, we defined the study area based on ecoregion boundaries to facilitate the identification of forest stability classes relative to a broad ecosystem type. In this way, we controlled for the existing differences in vegetation structure and floristics (e.g., dominant species, maximum canopy height and density) across boreal Russia, since these are largely influenced by the predominant physical environmental conditions. We used an ecoregion division for Russia from previous studies [51,52] that includes-from west to east-the administrative territories of Sverdlovsk, Tyumen, Omsk, Novosibirsk and Tomsk oblasts, Krasnoyarsk krai, and Irkutsk oblast. For the Amazon rainforest, the study region coincided with the Kayapó indigenous territory in the Brazilian states of Pará and Mato Grosso. This is an extensive area of largely primary tropical forest. We added a buffer zone of 60 km around the perimeter of the Kayapó territory to include areas that have been disturbed, degraded, and cleared by modern land use. The total study area comprised 120 million ha and 22 million ha of southern taiga and tropical forest, respectively.

Input Parameters
For each region, we used 16-year monthly time series (2003-2018) of (i) the fraction of photosynthetically active radiation intercepted by sunlit vegetation canopy (fPAR) and (ii) the shortwave infrared water stress index (SIWSI). Both metrics were obtained from MODIS satellite imagery at ~500 m (463 m) spatial resolution. Moreover, fPAR represents the fraction of incoming solar radiation between 400 and 700 nm absorbed by plant canopies and was taken from the MCD15A3H v6 Leaf Area Index/FPAR 4-day L4 global product [53]. SIWSI was derived from the near infrared (NIR; 841-876 nm) and the shortwave infrared (SWIR; 1628-1652 nm) bands, using the MOD09A1 v6 Terra Surface Reflectance

Input Parameters
For each region, we used 16-year monthly time series (2003-2018) of (i) the fraction of photosynthetically active radiation intercepted by sunlit vegetation canopy (fPAR) and (ii) the shortwave infrared water stress index (SIWSI). Both metrics were obtained from MODIS satellite imagery at~500 m (463 m) spatial resolution. Moreover, fPAR represents the fraction of incoming solar radiation between 400 and 700 nm absorbed by plant canopies and was taken from the MCD15A3H v6 Leaf Area Index/FPAR 4-day L4 global product [53]. SIWSI was derived from the near infrared (NIR; 841-876 nm) and the shortwave infrared (SWIR; 1628-1652 nm) bands, using the MOD09A1 v6 Terra Surface Reflectance 8-Day Global product, as follows [54]: For both products, we identified pixels that were affected by cloud contamination or sensor errors, using quality flags provided with the MODIS imagery. Particularly, we excluded pixels that contained high aerosol levels (bit 3), cirrus (bit 4), clouds (bit 5), or cloud shadows (bit 6). We then generated a mean monthly time series of fPAR and SIWSI from 2003 to 2018. We limited our analysis to June-August, which corresponds to the growing season in the southern taiga and the dry season in the Amazon rainforest. For the southern taiga, this procedure was useful to minimize variation related to seasonal changes in ground cover (i.e., snow) and vegetation greenness (i.e., phenology of deciduous trees). In contrast, the natural variation in satellite-derived vegetation indices over tropical rainforests is relatively small throughout the year due to the presence of a dense evergreen tree canopy. Therefore, inter-annual changes in forest greenness and canopy water content were expected to be particularly pronounced during the dry season (especially July and August) in the Amazon rainforest [55]. Finally, we masked out non-vegetative pixels corresponding to water, permanent wetlands, barren land, and urban areas, as well as croplands and crop-natural vegetation mosaics, using the MCD12Q1 Land Cover Type product corresponding to the year 2018 [56].
To assess forest-ecosystem stability, per-pixel metrics were generated from image stacks of mean June-August fPAR and SIWSI values over the period 2003-2018. For each variable, we first calculated the 16-year mean time-series value (fPAR mean and SIWSI mean ) and its coefficient of variation (fPAR CV and SIWSI CV ). These metrics provide information on the long-term average canopy greenness (fPAR), vegetation water stress (SIWSI), and their inter-annual variability (for mean and CV, respectively). The minimum seasonal fPAR (fPAR min ) and maximum seasonal SIWSI (SIWSI max ) values of the 16-year time series were calculated to provide information about the poorest (e.g., most water-stressed) vegetation condition attained by each pixel over the study period. In addition, we calculated linear annual trends (slopes) of seasonal fPAR and SIWSI for each pixel, using nonparametric Theil-Sen regression models. Slopes were obtained in order to characterize linear changes in vegetation status. Here, we focused on absolute slope values, which were indicative of the presence or absence of temporal trends but not of the direction of the trend (i.e., 'greening' or 'browning'). We considered a change in slope of both metrics as the direct result of forest disturbance due to either natural ecological processes, such as insect outbreaks, and/or human disturbances, such as logging and clearing. During the observation period, large positive slopes indicate forest succession following a prior disturbance, while large negative slopes indicate a disturbance during the time series. Areas with essentially no change were considered as those lacking alteration to surface features at the MODIS pixel scale (ca. 500 m).

k-Means Clustering of Pixels
Pixels with similar fPAR and SIWSI metric values were categorized into homogeneous classes by means of k-means clustering, which is an unsupervised algorithm that forms groups based on cluster centers. The use of an unsupervised algorithm has two main advantages compared to other approaches. First, it eliminates the need for extensive fieldbased observations to serve as training data for calibrating the remotely sensed metrics against surface features of interest. Second, this method is independent of particular environmental and forest characteristics, and, hence, it can be applied to any forested region, allowing for the application of a common approach to both boreal and tropical biomes. On the other hand, unsupervised algorithms require that great care is taken in selecting an appropriate area of interest and remotely sensed metrics to define the clusters since vague or inappropriate spectral grouping may not correspond to classes of interest. In addition, it demands a good knowledge of the ecosystems under evaluation to interpret unique classes produced during the classification, especially when they represent processes and not features on the ground [57].
All input variables underwent a two-step transformation prior to clustering. First, we tested whether the values of each variable were normally distributed through Shapiro-Wilk tests. If needed, the data were transformed to a new scale that better conformed to a normal distribution. The logarithmic transformation was used on all variables, except for slopes, where a cubic root transformation was applied. We then standardized the data to a uniform scale of 0-100, thereby allowing each variable to have the same weight in the clustering procedure. These procedures were performed at a biome scale for each region separately.
For the application of k-means, we defined 30 clusters (k = 30) that captured the natural variability at regional scales based on the previously defined parameters derived from fPAR and SIWSI time series, i.e., mean, minimum or maximum, CV, and absolute slopes for both fPAR and SIWSI. The k-value of 30 was chosen, as it provided a number of groups similar to the number of vegetation types described by the MCD12Q1 Land Cover Type product for both regions crossed with the three possible stability classes defined in our study (see details below). This number of groups was considered sufficient to distinguish between forest ecological conditions ranging from stable (including primary) forests to areas impacted by human land-use activities [11]. The clusters' centers were obtained separately for each region based on a sample of 10,000 random pixels (out of 6,332,993 pixels for the entire southern taiga ecoregion and 886,945 pixels for the Kayapó territory). Afterwards, every pixel within a region was assigned to its closest centroid by using Euclidean distances. The assignment was based on the highest probability of membership to a particular group. For each input variable, we subsequently obtained the mean values of the resultant clusters. Principal component analysis (PCA) was applied to the resulting multivariate matrix of mean values in order to summarize the variability captured by k-means clustering independently for each region.
Based on these mean values, the 30 groups were reclassified into the following three main categories, using a hierarchical clustering approach (Ward method) [58]: (i) forests with a high level of stability (HS), which would be representative of primary and mature forests; (ii) forests with a low level of stability (LS), which would include early succession and degraded stands; and (iii) non-forest (NF) areas, which would comprise deforested or naturally treeless patches. Since the initial 30 clusters were defined based on a sample of the total population, these categories represent areas (i.e., group of pixels) with the highest probability of membership to a given class. We refer to them as high-stability, low-stability, and non-forest areas for the sake of clarity. Overall, and based on previous evidence [55,59], we expected the following: • HS forests to be characterized by high fPAR and low SIWSI mean values over the time series and fPAR minimum and SIWSI maximum values close to the means, with high temporal stability of these metrics, i.e., low CV and slope; • LS forests to be characterized by lower fPAR and higher SIWSI mean values, more variable fPAR minimum and SIWSI maximum values, and less temporal stability, i.e., relatively high inter-annual fluctuation, high slope, or both; and • NF areas to display the lowest fPAR and highest SIWSI mean values, and also highly variable fPAR minimum and SIWSI maximum values, together with a large temporal instability, i.e., high CV and/or slope.
Spatial analyses, including image processing, parameter calculation, and k-means clustering, were performed by using Google Earth Engine (GEE) Javascript API [60]. Hierarchical (Ward) clustering was performed by using the 'stats' package in R software [61]. A workflow of the methodological approach to mapping forest stability is provided in Supplementary Figure S1a.

Validation Framework
We validated our forest ecosystem stability spatial data layer for each study region with three independent sources of data and approaches. It is important to note, however, that currently available data products did not allow us to directly validate the low-stability (LS) class, precisely because our method is the first to categorize LS forests at landscape scales, using remote-sensing time series. First, the accuracy of our product was visually assessed against time series of remotely sensed reference imagery at randomly sampled locations within the study regions [62]. The reference dataset consisted of Landsat TM/ETM+/OLI imagery, which was used by an independent expert to establish whether a forest area within a MODIS pixel footprint had been disturbed during the past several decades [63]. This analysis covered the period from 1984 to 2018 and enabled tracking forest disturbances before (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) and during the study period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). Poor pixels with low image quality or high cloud and shadow contaminations, which resulted in an insufficient level of certainty given by the expert (<50%), were removed from the analysis. The validation analysis focused on detecting changes in forest structure (i.e., tree canopy cover) to differentiate stable forests from forests with altered conditions. Areas of stable or disturbed forests were validated by using a probability-based stratified random sampling per biome. Our forest stability maps were used to create two major sub-strata per biome: no change, which corresponded to high stability forests (HS); and change, which was defined as highly disturbed, deforested, or non-forest (NF) areas. Low-stability (LS) forests were excluded from this analysis since this class represents more subtle changes in remotely sensed metrics of fPAR and SIWSI and often do not correspond to features that are visually detectable by using optical imagery (hence emphasizing the benefit of our approach). In total, 545 and 320 pixels were evaluated for the Southern Siberian taiga and the Kayapó territory. Image interpretation of time-series imagery was performed by using the Google Earth Timelapse app (https://earthengine.google.com/timelapse/; accessed on 19 November 2021).
Next, we compared stability classes with stand-disturbance history retrieved from a number of independent datasets over the period 2003-2018, as described below. For the southern taiga, we manually digitized areas with substantial forest-cover change corresponding to two major disturbances (fire and logging) over the period 2003-2018, using high-resolution Landsat imagery (ca. 30 m spatial resolution). The geographic area for validation comprised the Angara region, a deforestation hot spot in Siberia resulting from timber harvest and recurrent fires [50] that covers 12.7 million ha (i.e.,~10% of the study region). We used 30 m resolution data in order to retrieve information at a MODIS sub-pixel level, with the purpose of characterizing the relatively small areas of cut blocks and the heterogeneity of fire impacts. Fractions of MODIS pixels affected by either timber harvest or fire were generated separately. Thresholds of 10% and 90% were applied to these pixel fractions in order to identify areas with varying disturbance intensities: (i) undisturbed areas (<10% of each pixel affected by fire or logging), (ii) moderately disturbed areas (between 10% and 90%), and (iii) highly disturbed areas (>90%).
For the Kayapó territory, we used a recently updated deforestation spatial layer from the Brazilian space agency (PRODES and INPE), which is based on analysis of the Landsat archive (30 m spatial resolution) and provides information on the distribution of primary forests, degraded forests (mainly due to logging) and non-forest areas [64]. Using these data, we estimated the proportion of each of these categories within each MODIS pixel.
Finally, we evaluated the extent to which the mapped forest stability classes in both regions were consistent with estimates of tree canopy cover change, as reported in Reference [65]. This dataset provides information on areas experiencing tree cover loss or gain globally at ca. 30 m resolution. We evaluated the proportion of cumulative forest change over the period 2003-2018 within each MODIS pixel across the study regions. We applied the following threshold values: no change in tree cover (<10% of each pixel showing tree cover change), moderate change (between 10% and 90%), and substantially altered tree-cover areas (>90%). A summary of the datasets used for validation of the final forest stability maps is provided in Supplementary Figure S1b.

Results
The hierarchical clustering approach identified three classes based on MODIS time series of fPAR and SIWSI metrics for the southern taiga in Siberia (Figure 2a) and the Amazon rainforest in the Kayapó territory (Figure 3a). These classes captured 65% and 71% of the total variability present in the 30 k-means groups for Siberia and Kayapó, respectively. A PCA performed on the 30 k-means groups returned two principal components (PCs), which accounted for 86% and 88% of the total variance for the Southern Siberian taiga (Figure 2b) and the Kayapó territory (Figure 3b), respectively. These analyses unequivocally led to the definition of distinct classes of forest stability in both regions that are well illustrated by the loadings of each classificatory variable in the biplots shown in Figures 2b and 3b. Particularly, high-stability (HS) forests were characterized by conditions representative of productive and stable forest ecosystems: high fPAR (mean and min), low SIWSI (mean and max), low CVs, and low absolute slopes (Figures 2 and 3; Table 1). In contrast, low-stability (LS) forests were characterized by conditions of lower productivity and stability, including lower fPAR (mean and min), higher SIWSI (mean and max), higher CVs, and higher slope magnitudes (moderately higher in Siberia and much higher in Brazil; Figures 2 and 3; Table 1). Finally, non-forested (NF) areas displayed the lowest fPAR (mean and min), the highest SIWSI (mean and max), and were much less stable than HS (higher CVs and slope magnitudes; Figures 2 and 3; Table 1 biplots shown in Figures 2b and 3b. Particularly, high-stability (HS) forests were characterized by conditions representative of productive and stable forest ecosystems: high fPAR (mean and min), low SIWSI (mean and max), low CVs, and low absolute slopes (Figures 2 and 3; Table 1). In contrast, low-stability (LS) forests were characterized by conditions of lower productivity and stability, including lower fPAR (mean and min), higher SIWSI (mean and max), higher CVs, and higher slope magnitudes (moderately higher in Siberia and much higher in Brazil; Figures 2 and 3; Table 1). Finally, non-forested (NF) areas displayed the lowest fPAR (mean and min), the highest SIWSI (mean and max), and were much less stable than HS (higher CVs and slope magnitudes; Figures 2 and 3; Table  1   The resulting regional forest stability layers revealed distinct spatial patterns that were consistent with other geospatial products, including higher-resolution imagery, digitized disturbance data, regional and global land-cover maps, and layers of tree-cover changes (see Section 2.3 for details). In particular, a visual comparison of forest stability maps with higher-resolution satellite images revealed a high level of agreement between stable forest canopies and areas clearly affected by human impacts (Figure 4c,d). The results of the visual expert assessment of structural changes in tree canopy cover over the period 1984-2018 showed a high accuracy of characterization of forest stability classes in both biomes ( Table 2). In the southern taiga, the lowest user's accuracy of 76%, reflecting commission error, was found in NF areas, and the lowest producer's accuracy of 91%, reflecting omission error, was found for HS forests. In the Kayapó territory, both user's and producer's accuracies were above 97%. The overall accuracy was 92% for the boreal and 99% for the tropical biomes (Table 2). It is worth noting that our approach was able to     commission error, was found in NF areas, and the lowest producer's accuracy of 91%, reflecting omission error, was found for HS forests. In the Kayapó territory, both user's and producer's accuracies were above 97%. The overall accuracy was 92% for the boreal and 99% for the tropical biomes (Table 2). It is worth noting that our approach was able to capture the long-term effect of past disturbances occurring up to two decades before the study period (i.e., 1984-2002 period).    Table 2. Estimated error matrix for expert validation of the resultant forest stability maps for southern taiga ecoregion (Siberia, Russia) and Kayapó territory (Amazon basin, Brazil). The validation is based on expert knowledge of the forested landscapes through the visual interpretation of changes in forest cover by using high-resolution imagery over the period 1984-2018. Forest stability classes are in rows, while the reference classes are in columns. This visual assessment was performed by using the Google Earth Timelapse app. All pixels for which expert certainty was low (<50%) were removed from the analysis. HS = high-stability forests; LS = low-stability forests; NF = non-forest areas.

Reference Data
User's Accuracy (%) Importantly, a pixel-by-pixel comparison of our stability map to ten alternative runs of our approach by using different seeds for the random pixels used to train the k-means algorithm yielded consistent outputs (Supplementary Figure S3). Particularly, the lowest probability of a pixel to be assigned to a particular stability class was above 60% in both regions. Overall, the degree of agreement between our stability map and the alternative maps was 93% and 98% in the southern taiga and the Kayapó territory, respectively. Among the three stability classes, the highest degree of agreement was found for the high-stability forests (HS; 96% in the southern taiga and 99% in the Kayapó territory), followed by non-forest (NF) area (89% in the southern taiga and 100% in the Kayapó territory) and low-stability (LS) forest (90% in both regions).
Validation data on known fire and logging disturbances further confirmed the accuracy of the approach for the Angara region in Siberia ( Figure 5). The majority of pixels classified as HS corresponded to a total lack of disturbance by either fire (79%) or logging (92%) during the study period. Pixels with intermediate levels of disturbance (up to 90%) corresponded mainly to LS (50%), and pixels showing high disturbance (>90%) were mainly represented by NF (56%). In the Kayapó territory, the vast majority of area mapped as primary forests according to the PRODES deforestation data was classified as HS ( Figure 6), whereas areas mapped as deforested, non-forest, and other were dominated by LS and NF stability classes.
Among the three stability classes, the highest degree of agreement was found for the highstability forests (HS; 96% in the southern taiga and 99% in the Kayapó territory), followed by non-forest (NF) area (89% in the southern taiga and 100% in the Kayapó territory) and low-stability (LS) forest (90% in both regions).
Validation data on known fire and logging disturbances further confirmed the accuracy of the approach for the Angara region in Siberia ( Figure 5). The majority of pixels classified as HS corresponded to a total lack of disturbance by either fire (79%) or logging (92%) during the study period. Pixels with intermediate levels of disturbance (up to 90%) corresponded mainly to LS (50%), and pixels showing high disturbance (>90%) were mainly represented by NF (56%). In the Kayapó territory, the vast majority of area mapped as primary forests according to the PRODES deforestation data was classified as HS (Figure 6), whereas areas mapped as deforested, non-forest, and other were dominated by LS and NF stability classes.  For additional validation, we compared our forest stability classes to a widely used dataset of global changes in tree canopy cover [65]. Consistent with our other comparisons, a majority of areas mapped as low forest loss or gain were classified as HS (58% and 72% for the southern taiga and Kayapó territory, respectively), whereas areas mapped as substantial forest loss or gain (change detected for >50% of a single MODIS pixel) were classified primarily as LS (30% and 60%) or NF (54% and 40%) (Supplementary Figure S4).
Overall, the HS class corresponded to forest patches having a relatively high tree canopy cover with minimal disturbance, which was the most prevalent category in Southern Siberia (54% of total land; Figure 4a) and in the Kayapó territory (63% of total land; Figure 4b). LS encompassed both degraded forests impacted by fires and human activity (e.g., logging, mining, or agricultural practices) over the study period ( Figures 5 and 6 and Supplementary Figure S3), as well as areas with naturally sparse vegetation cover. The latter was the case for the rear (i.e., southernmost) edge of the Southern Siberian taiga, which was characterized by a relatively low tree-cover density in the ecotone between the southernmost taiga and the forest steppe (Supplementary Figure S5a,c). LS accounted for 30% of the area in Southern Siberia and 19% in Kayapó. NF targeted areas with very low treecover density (Figure 1). These included natural treeless patches (e.g., grasslands or rocky areas) or, alternatively, deforested zones resulting from land conversion to agricultural fields, clear-cut or severe fires (Figures 5 and 6 and Supplementary Figures S3 and S4), resulting in significant temporal (e.g., fires or logging) or permanent losses of tree cover (e.g., infrastructure development or agriculture). NF pixels accounted for 16% of the total area in Southern Siberia (Figure 4a) and 18% in Kayapó (Figure 4b). For additional validation, we compared our forest stability classes to a widely used dataset of global changes in tree canopy cover [65]. Consistent with our other comparisons, a majority of areas mapped as low forest loss or gain were classified as HS (58% and 72% for the southern taiga and Kayapó territory, respectively), whereas areas mapped as substantial forest loss or gain (change detected for >50% of a single MODIS pixel) were classified primarily as LS (30% and 60%) or NF (54% and 40%) (Supplementary Figure S4).
Overall, the HS class corresponded to forest patches having a relatively high tree canopy cover with minimal disturbance, which was the most prevalent category in Southern Siberia (54% of total land; Figure 4a) and in the Kayapó territory (63% of total land; Figure 4b). LS encompassed both degraded forests impacted by fires and human activity (e.g., logging, mining, or agricultural practices) over the study period ( Figures 5 and  6 and Supplementary Figure S3), as well as areas with naturally sparse vegetation cover. The latter was the case for the rear (i.e., southernmost) edge of the Southern Siberian taiga, which was characterized by a relatively low tree-cover density in the ecotone between the southernmost taiga and the forest steppe (Supplementary Figure S5a,c). LS accounted for 30% of the area in Southern Siberia and 19% in Kayapó. NF targeted areas with very low tree-cover density (Figure 1). These included natural treeless patches (e.g., grasslands or rocky areas) or, alternatively, deforested zones resulting from land conversion to agricultural fields, clear-cut or severe fires ( Figures 5 and 6 and Supplementary Figures S3 and  S4), resulting in significant temporal (e.g., fires or logging) or permanent losses of tree cover (e.g., infrastructure development or agriculture). NF pixels accounted for 16% of the total area in Southern Siberia (Figure 4a) and 18% in Kayapó (Figure 4b).

Discussion
Here we defined forest stability classes derived from remote-sensing time-series data related to photosynthetic activity (fPAR) and water stress (SIWSI). These metrics are closely associated with ecophysiological processes that are indicative of photosynthetic

Discussion
Here we defined forest stability classes derived from remote-sensing time-series data related to photosynthetic activity (fPAR) and water stress (SIWSI). These metrics are closely associated with ecophysiological processes that are indicative of photosynthetic activity and vegetation structure at the forest-stand and landscape levels [55,59], and the source data are readily available globally. The selection and interpretation of the metrics are informed by ecological theories of resource optimization, microclimate buffering, biological redundancy, and adaptive capacity [66,67]. We focused on representative regions of boreal and tropical forest because the majority of remaining primary forests are located in these biomes [68,69].
A large-scale validation dataset against which to measure the accuracy of our stability maps is not yet available, and particularly based on long-term ground truth data. Thus, to evaluate the accuracy of our approach, a few comparisons were made with other existing regional or global datasets which employed high-resolution data sources. Although these cannot fully substitute for validation data, they help to characterize the maps by yielding a degree of agreement between products that were derived entirely independently from one another. The results of the validation analyses suggest that our approach is applicable and sufficiently accurate to map forest stability (in terms of constancy) within boreal and tropical forest ecoregions. Both spatial and temporal assessments of forest stability classes showed consistent performance of our approach. The distinction between high-stability forests and non-forest areas showed a high level of agreement with other sources. On the other hand, the low stability class is, by its nature, difficult if not impossible to independently verify using existing products and optical imagery. This forest class is not necessarily subject to large-scale or intense disturbance but is, instead, more variant in time and susceptible to changing environmental conditions. Such dynamics can only be detected using long-term time series, a benefit of our approach compared to existing methods focused on distinct periods in time. It should also be noted that our digitized fire polygons in Siberia may represent either stand-replacing or non-stand-replacing fires, and that future work could further explore the impact of different fire severity levels on forest stability.
Our approach complements other approaches to assessing forest conditions globally. These include mapping intact forest landscapes [41], human pressure indicators [70][71][72][73], habitat fragmentation [74], tropical hinterland forests [22,75], and vegetation structural condition [76,77]. Each of these mapping techniques has its own advantages and limitations. For example, existing intact forest landscape layers provide information on globally consistent forested landscapes with little to no human interference [41], but they are intended to identify large forest masses (i.e., the associated patch size thresholds exclude valuable primary forests smaller than 50,000 ha), and the identified patches include non-forest land cover. Global evaluations of forest loss using standardized definitions can be used at a pixel level and can be useful for a wide variety of applications [65]. However, such definitions, based on tree canopy cover without complementary analysis of forest ecological conditions, equate heavily managed forests and tree plantations with primary forests, leading to potential misinterpretations and misuse [78].
Our approach to mapping forest stability bypasses some of these limitations by not relying on proxies of human influence, not requiring intensive machine-learning parameterization and training, and integrating multi-year information on time-series dynamics. It is also applicable across forest biomes and allows for the mapping of both large tracts and small patches of stable forests at any scale down to a minimum area consistent with that found necessary to maintain forest interior micro-environmental conditions (generally 500-2000 m) [79,80]. Smaller areas of stable forest constitute a substantial proportion of the natural forest cover in many regions [81], which are of considerable conservation value, including refugia habitat, biodiversity, and source populations for landscape restoration, among other ecosystem services [82][83][84][85][86].
While MODIS satellite imagery has coarser spatial resolution as compared with widely used Landsat data, we argue that this scale is appropriate for a landscape-level analysis of forest stability. By mapping the constancy property at a 500 m pixel resolution (~25 ha), we capture gradients in stability caused by edge effects, which have been estimated at over 1000 m for wet tropical forests (suggesting an absolute minimum area for a patch of undisturbed forest to be around 100 ha) [13,55] and 100-400 m for boreal and temperate forests (and a minimum patch of undisturbed forest being 1-16 ha) [87][88][89]. Boreal forest landscapes are naturally characterized by patches of different forest ecosystem types, including multiple seral stages following wildfire [90] which have been found to range in size from 2 to 650 ha. The spatial resolution of our analysis, therefore, is commensurate with the scaling of forests with a high level of stability in tropical biomes and for larger sized patches in boreal biomes. Likewise, by leveraging multi-year time series, we were able to capture the temporal dimension of forest degradation, which can take several years to manifest [91].
Our approach is limited in a number of important aspects that require further attention. The stability metric used was designed to address only the constancy property of stability, defined here as no major changes in forest canopy condition over a 16-year period (2003-2018). To evaluate the other two stability properties of resilience and persistence [30], consideration needs to be given to analyzing multi-patch dynamics and landscape-level metrics over a longer time series. For example, the available data are unable to represent resilience following state-altering disturbances and the potential for ecosystems to return to their pre-disturbed state or the optimum operating point given prevailing conditions (see Reference [92], for example). This is particularly relevant to boreal forests, where recurrent wildfires are part of their natural disturbance regime, and resilience needs to be considered from a landscape-wide perspective to take into account the shifting mosaic of seral stages over time. This is in contrast to tropical forests, which are characterized by a continuous closed canopy (i.e., projected foliage cover >70%) [93]. Thus, an expanse of structurally intact boreal forest canopy is unequivocally a stable forest, but patches of disturbed seral stages may be the result of either natural disturbances, especially wildfire, or human landuse impacts. On the other hand, loss of a structurally intact canopy in the wet tropics is a more reliable indicator of human rather than natural disturbance (notwithstanding impacts from large blowdowns [94] and tropical storms [95]).
Our stability metric did not distinguish between natural and anthropogenic disturbances. These contrasting dynamics could potentially be characterized by broken-line (segmented) regression methods [96], especially if progressively long time series of remotely sensed data are available. Uptrends and downtrends can also be used to help differentiate between vegetation greening and browning. In the tropics, longer series may particularly aid in disturbance detection since wet forests can reach a closed tree canopy within a few decades after disturbance. In the boreal, however, a second stage of analysis using finer spatial resolution satellite data and object-based image classification techniques [97] could potentially tease out human land-use impacts based upon their characteristic spatial patterning compared to those arising from natural disturbances [98]. Therefore, our approach is a first step towards mapping all three stability properties at a landscape scale [99].
Within a given ecoregion, our method also tends to classify naturally low-canopycover ecosystems as low stability compared with neighboring, relatively closed-canopy forests. The high-stability class is oriented, therefore, towards the more closed-canopy forests at an ecoregion level. Although the total area of lower-canopy ecosystems is minor in the southern taiga, possible variants of this methodology could employ metrics of habitat quality tailored to specific ecological units (e.g., land cover types), decouple fPAR magnitudes from stability metrics, or be mapped at a higher spatial resolution relevant to local-scale assessments, including through fusing of MODIS data with higher-spatialresolution imagery, such as Landsat and Sentinel-2 (see Reference [100], for example). Finally, we note that this work has focused on two case-study forests that are typical of boreal and wet tropical forest ecosystems. Thus, further testing of the approach is needed in temperate and seasonally dry tropical forests.
Information on stable forests-including the constancy metric applied here-is relevant to national and international policies concerning biodiversity conservation and climate-change mitigation. At a national level, information about high-stability forests can be used to inform systematic conservation planning [101]. Information on the location of less stable forests is also relevant for ecological restoration, reforestation, and proforestation activities [102][103][104][105][106], or forestry intensification to alleviate logging pressure on primary forests [7,8,105]. Moreover, from a spatial ecology perspective, patches of disturbed forests are an inherent component of most temperate and boreal forest landscapes (see Reference [107], for example) and are valuable for biodiversity and the persistence of particular species dependent on early successional communities. In general, as long as enough of the surrounding forests are stable (i.e., resilient to disturbances-be they natural or anthropogenic), such ecosystems are sustainable [108].

Conclusions
Our methodology represents an ecophysiologically based approach to monitor and map ecosystem stability, as showcased for two contrasting forest biomes with the potential to be applied elsewhere. If upscaled and implemented globally, it would provide a uniform, objective source of information on a key property of forest ecosystems that is consistent with widely accepted definitions of primary forests. Recognizing the limitations of the approach, the data generated could be used as an independent source to cross-validate other local or global forest assessments. Because it is designed to be repeatable based on objective statistics with minimal expert intervention, assessments of global forest stability could be regularly updated by extending the time series as new satellite data become available. The information could be used to inform national and international policy settings that aim to prioritize conservation investments targeting stable forest ecosystems to maximize outcomes for biodiversity conservation, carbon storage, and other ecosystem services. Ecoregions that are still relatively free from the damaging impacts of large-scale human activities are rapidly disappearing, yet they provide valuable ecosystem service benefits to people. Because effective restoration can be costly and requires long-term management, identifying actions that can retain stable forests, including primary forests and intact forest landscapes, is a recommended strategy for advancing international commitments for conservation, biodiversity, and climate mitigation. Our approach, which maps the ecosystem stability of forests in terms of their constancy, contributes to the growing information based on the condition of forest ecosystems that remains an underassessed component of current global data layers.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/rs14153813/s1. Figure S1: Methodology flowchart and validation procedure for mapping and accuracy assessment of forest ecosystem stability within major biomes. Figure S2: Representative time series of mean fPAR and SIWSI for southern taiga ecoregion (Siberia, Russia) and Kayapó territory (Amazon basin, Brazil) over the period 2003-2018. Figure S3: Spatial constancy analysis for southern taiga ecoregion (Siberia, Russia) and Kayapó territory (Amazon basin, Brazil). Figure S4: Spatial distribution of forest cover changes for southern taiga ecoregion (Siberia, Russia) and Kayapó territory (Amazon basin, Brazil). Figure S5: Spatial distribution of vegetation types from the Land Cover Type product corresponding to the year 2018 for southern taiga ecoregion (Siberia, Russia) and Kayapó territory (Amazon basin, Brazil). Table S1: Major climate characteristics of southern taiga ecoregion (Siberia, Russia) and Kayapó territory (Amazon basin, Brazil).