Spatial patterns of logging-related disturbance events: a multi-scale analysis on forest management units located in the Brazilian Amazon

Selective logging has been commonly mapped using binary maps, representing logged and unlogged forests. However, binary maps may fall short regarding the optimum representation of this type of disturbance, as tree harvest in tropical forests can be highly heterogeneous. The objective of this study is to map forest disturbance intensities in areas of selective logging located in the Brazilian Amazon. Selective logging activities were mapped in ten forest management units using Sentinel-2 data at 10 m resolution. A spatial pattern analysis was applied to the logging map, using a moving window approach with different window sizes. Two landscape metrics were used to derive a forest disturbance intensity map. This map was then compared with actual disturbances using field data and a post-harvest forest recovery analysis. Disturbed areas were grouped into five distinct disturbance intensity classes, from very low to very high. Classes high and very high were found to be related to log landings and large felling gaps, while very low intensities were mainly related to isolated disturbance types. The post-harvest forest recovery analysis showed that the five classes can be clearly distinguished from one another, with the clearest differences in the year of logging and one year after it. The approach described represents an important step towards a better mapping of selectively logged areas, when compared to the use of binary maps. The disturbance intensity classes could be used as indicators for forest monitoring as well as for further evaluation of areas under forest management.


Introduction
Selective logging is a major anthropogenic activity and a key component of the forest degradation process in the tropics. In the Brazilian Amazon, selective logging activities have been reported on areas as large as those reported as deforested . The impacts of selective logging have been documented for a wide variety of biotic and abiotic indicators (Meijaard et al. 2005;Olander et al. 2005;Asner et al. 2009;Burivalova et al. 2014;Darrigo et al. 2016;de Carvalho et al. 2017;Stas et al. 2020). While there are clear signs of negative impacts associated with selective harvest operations, this activity is less deleterious to the environment than deforestation, which is a complete land cover and land use change. In fact, it has been argued that selectively logged forests, if properly managed, can maintain important biodiversity values, carbon stocks and other ecosystem services, unlike deforested areas (Putz et al. 2012;Laurance and Edwards 2014;Edwards et al. 2019). Therefore, taking into account that selective logging is the main harvest technique employed in natural tropical forests (Poudyal et al. 2018), efforts to ensure its sustainability are of great interest.
In this context, the ''sustainable management of forests'' was proposed to be included in the United Nations Framework Convention on Climate Change (UNFCCC) Reducing Emissions from Deforestation and Forest Degradation (REDD?) initiative in 2007, and was finally approved in 2013 (Seymour and Busch 2016). At the same time, several countries have made substantial progress towards the implementation of forest management activities in the tropics (Blaser et al. 2011). While there is still a debate about the use of the word ''sustainability'' (Putz 2018;Tegegne et al. 2018), it has been reported that forest management in tropical forests may help to stall deforestation activities (Burivalova et al. 2020). Concurrently, in recent decades policies such as the Forest Law Enforcement, Governance and Trade (FLEGT) policy of the European Union were created to foster the sustainability and legality of tropical forest products (Tegegne et al. 2018). The correct implementation of forest management activities, ensuring their sustainability, and the legality of forest products, are highly dependent on the countries' monitoring capacity. In fact, it has been reported that the absence of monitoring and external control favors unsustainable practices (Rutishauser and Herold 2017). This highlights the importance of reliable monitoring systems as tools to promote the better use of forests.
Efforts to develop Criteria and Indicators (C&I) for monitoring sustainable forest management (SFM) activities in tropical forests have been made by a variety of national and international organizations (Elias 2004;ITTO 2016;Linser et al. 2018). Many indicators rely on field data; others, however, can only be measured with remote sensing techniques, due to the fact that many SFM plans implemented in tropical forests, especially in the Amazon basin, are located in areas that are difficult to access. Therefore, information derived from satellite images are used as a proxy to measure the impacts of selective logging over time and covering large areas. This highlights the importance of investigating the use of different satellite input data in landscape ecology studies.
In the last two decades, we have seen a major development of remote sensing techniques applied to the assessment of selective logging in the Brazilian Amazon Matricardi et al. 2010;Souza Jr. et al. 2013;Pinagé et al. 2016;Pinheiro et al. 2016;Tritsch et al. 2016;Grecchi et al. 2017;Dalagnol et al. 2019;Hethcoat et al. 2019;Lima et al. 2019;Shimabukuro et al. 2019;Bullock et al. 2020). While the majority of these studies focused on the mapping of selective logging, some of them assessed the intensity of forest disturbances and/or their spatial pattern either using a grid-based approach (Pinheiro et al. 2016;Grecchi et al. 2017) or using land tenure databases (Tritsch et al. 2016). However, there is a need to investigate patterns of forest disturbances caused by selective logging in forest management units (FMU), rather than in artificial landscape units such as regular grids. This is the case because the FMU is a unit where forest management activities are actually happening and, therefore, it is expected that most of the forest disturbances within these areas are caused by selective logging activities. While within a grid, arbitrarily placed in the landscape, many types of forest disturbances could be mapped. An FMU can be interpreted as ''a clearly defined forest area managed to a set of explicit objectives according to a long-term management plan'' (ITTO 2016). Investigating actual disturbance patterns in government authorized FMUs is crucial for understanding how the SFM has been implemented. The few existing studies assessing forest disturbances within FMUs cover only forest concessions (Pinagé et al. 2016;Dalagnol et al. 2019;Hethcoat et al. 2019), but not smaller, private forest properties. Moreover, none of these studies had the objective of mapping forest disturbance intensities. Therefore, a detailed spatial analysis of forest disturbances caused by selective logging is necessary.
So far, the vast majority of remote sensing studies have focused on mapping the extent of selectively logged forests as binary maps of disturbed/undisturbed forests Matricardi et al. 2010;Souza Jr. et al. 2013;Pinheiro et al. 2016;Tritsch et al. 2016;Grecchi et al. 2017;Lima et al. 2019;Shimabukuro et al. 2019). Negative impacts associated with selective logging in the tropics can include the direct loss of biomass, damage to the remaining trees, compromising the natural regeneration, soil compaction and susceptibility of the remaining forest to fires, due to the changes in the microclimate associated with the canopy openings (Meijaard et al. 2005;Asner et al. 2009;Poudyal et al. 2018). However, these impacts vary according to the type of logging infrastructure created during the harvest process. Secondary roads, log landings, felling gaps and skid trails may impact forests in a different manner, with some of these features, such as log landings and logging roads, causing much more impact to the landscape than the others (Asner et al. 2002;Pinagé et al. 2016;de Carvalho et al. 2017). Therefore, binary maps measuring the extent of selectively logged areas as an indicator of forest disturbances may not be sufficient to describe the logging impacts in areas under SFM. A spatial pattern analysis aiming at the mapping of forest disturbance intensities could address this research gap.
The main premise of this study is that selective logging activities in tropical forests are highly heterogeneous and that different types of patches, with distinct disturbance intensities, are produced during harvest operations. A patch implies a ''relatively discrete spatial pattern'' and a ''relationship of one patch to another in space'', considering the surrounding affected and unaffected areas (White and Pickett 1985). Therefore, the type of spatial configuration of forest disturbance patches formed during selective logging processes could give us important information about the selective logging impact in a given area. This information is crucial for practical reasons, e.g. the future improvement of C&Is for SFM. A spatial pattern analysis (Riitters 2019) can be used to create a map accounting for the intrinsic heterogeneity within the areas of selective logging. Such spatial pattern analysis allows categorizing each pixel of a binary map into discrete disturbance intensity classes while keeping the spatial resolution of the pixel Riitters et al. 2017).
In the present study, we analyzed the spatial patterns of forest disturbances caused by selective logging within forest management units located in the Brazilian Amazon. The main objective is to derive a map with forest disturbance intensity classes, based on the two fundamental landscape metrics: composition and configuration. In particular, we addressed the following research questions: (i) How do the values of composition and configuration change over multiple observation scales? (ii) What would be the observation scale that could best distinguish between different types of disturbances associated with selective logging? (iii) Is it possible to associate forest disturbance intensity classes, derived from spatial pattern analysis, with actual disturbances caused by selective logging? (iv) How different is the post-harvest recovery process among distinct classes of forest disturbance intensity?

Methods
The methodology is made up of three main steps: first, logging activities were mapped using high resolution Sentinel-2 imagery (Drusch et al. 2012). The result is a binary map containing two classes: disturbed and undisturbed forests. Second, a spatial pattern analysis was applied to this binary map. Third, landscape metrics were used to derive a forest disturbance intensity map via a cluster analysis. The method used to generate the final forest disturbance intensity map is partially based on the multi-scale approach proposed by Zurlini et al. (2006) and refined by Riitters et al. (2017). In this study the multi-scale approach was applied in a particular data set with a finer spatial resolution (10 m), making it particularly suitable to capture the heterogeneity present in areas of selective logging. By adapting this method to small-scale logging-related disturbances, it was possible to determine five categories of logging and, therefore, go beyond binary maps.

Study area and selection of the forest management units
The study area is located in the southern region of the Brazilian State of Amazonas. Amazonas is the largest state in Brazil (1.5 million km 2 ) and holds the largest area of intact forests within the Brazilian Amazon (INPE 2020). The region has a tropical monsoon climate (Am) according to the Köppen climate classification system (Alvares et al. 2013). Considering the TerraClimate data set (Abatzoglou et al. 2018) for the last 40 years , the averaged minimum and maximum air temperature are 21.0°C and 31.5°C, respectively, and the annual rainfall ranges from 2219 to 3370 mm. The short dry season (precipitation less than 60 mm per month) occurs between June and August (Abatzoglou et al. 2018). The study area is covered predominantly with terrafirme forests (non-flooded) (IBGE 2012). The study was carried out in a focus area near the Transamazon Highway (BR 230), in the district of Santo Antônio do Matupi (Fig. 1). This is currently one of the main timber production zones within the State of Amazonas (IPAAM 2020), while the surrounding area is a deforestation hotspot (INPE 2020).
The study was carried out in FMUs licensed by the Institute of Environmental Protection of the Amazonas State (IPAAM) under the category of ''sustainable forest management of major impact'' (CEMAAM 2013). Therefore, hereafter, the term ''sustainable forest management -SFM'' follows its definition under Brazilian legislation (CONAMA 2009;CEMAAM 2013). In the State of Amazonas, logging intensity is limited to 25 m 3 /ha, within areas suitable for timber extraction within each SFM area (CEMAAM 2013). However, in general, the logging intensity is much lower than 25 m 3 / ha (Lima et al. 2019).
According to data obtained from IPAAM, a total of 130 SFM areas received authorization for logging between 2016 and 2017, considering all municipalities in the State of Amazonas. From these areas, 24 (18.4% of the licensed areas) are located within the Area of Interest (AOI). Of these, 10 areas that were selectively logged in 2017 were chosen for the spatial pattern analysis. Therefore, in summary, these areas were selected based on their location (surrounding the Village of Santo Antônio do Matupi) and also because they had experienced harvest activities during the dry season of 2017, the target year of analysis. Moreover, field data regarding logging activities were collected Fig. 1 Location of the study area and the forest management units (black polygons, numbered from 1 to 10). The background map with land cover classes (year 2017) was retrieved from PRODES (INPE 2020) in seven FMUs (Lima et al. 2019), which allowed us to make further analysis and comparisons.

Spatial pattern analysis
The spatial pattern analysis was carried out using a binary map of disturbed and undisturbed forests produced for the period 2016-2017. This map was derived using Sentinel-2 images at 10 m resolution. It was based on a change detection approach (Langner et al. 2018), which was already tested for Sentinel-2 data in the some of the FMUs investigated in the current study (Lima et al. 2019). Sentinel-2 is a suitable satellite to map small-scale forest disturbances caused by selective logging (Lima et al. 2019). The steps to produce this binary map are described in the Supplementary Material. A variety of metrics exists to estimate spatial patterns of forest disturbances (Turner and Gardner 2015). However, ultimately, the fundamental elements of landscape patterns are described by composition and configuration (Li and Reynolds 1994;Gustafson 2019;Riitters 2019).
Composition can be calculated as the proportion of the map occupied by the land cover class of interest (Pd = proportion). Therefore, pixels belonging to the disturbance class can have Pd values ranging from 0 to 100% (Vogt 2019a). An example of a Pd calculation for a square moving window of 7 9 7 pixels is shown in Fig. 2. Configuration can be described by a landscape metric known as adjacency (Riitters 2019). Adjacency, also called contagion (Riitters et al. 1996), is used to distinguish between landscape patterns that are clumped or dispersed (Turner and Gardner 2015). Here, adjacency is described as the conditional probability (using the 4-neighbor rule) that a focal class pixel is adjacent to another focal class pixel (Pdd = adjacency). Pixels belonging to the disturbance class are also going to have Pdd values ranging from 0 to 100% (Fig. 3). Here, adjacency is defined as a ''class-level contagion'' because just one land cover class, the disturbance class, was used to calculate Pdd values (Vogt 2019b). Both metrics were calculated using a square moving window with a variety of sizes, using ''GuidosToolbox'', a software developed by the European Commission (Vogt and Riitters 2017).
The spatial variation of a system property will depend on the spatial scale at which the property is measured and the size of the mapping unit (Gustafson 1998). In this study, the size of the mapping unit (also called the grain size) is the 10 m spatial resolution of the input map, derived from Sentinel-2 data. The spatial scale follows the conceptual model proposed by Zurlini et al. (2006), where the scale of observation is defined by a fixed-area following square moving windows of different sizes. Therefore, Pd and Pdd were calculated for 8 fixed-area windows: 3 9 3 pixels (0.09 ha), 5 9 5 (0.25 ha), 7 9 7 (0.49 ha), 9 9 9 (0.81 ha), 11 9 11 (1.21 ha), 15 9 15 (2.25 ha), 21 9 21 (4.41 ha) and 25 9 25 (6.25 ha). Hereafter, the term ''landscape extent'' is used to Fig. 2 Example of the proportion (Pd) metric computation from a binary map. Left: Study area No. 6, binary map of undisturbed/ disturbed forests; Center: Moving window (7 9 7 pixels) surrounding the focal pixel (in blue); Right: Proportion (Pd) map, containing pixels with the proportion of the surrounding landscape encompassing disturbed pixels. In the case of the focal pixel exemplified Pd = 76% define different moving window sizes and when the fixed-area windows are described as 3 9 3, for example, it means that the size of the moving window is 3 9 3 pixels. Therefore, a moving window of 3 9 3 pixels can also be called a landscape extent of 3 9 3 pixels. Previous studies have used larger window sizes to describe landscape patterns at multiple scales Riitters et al. 2017;Wardlaw et al. 2018). However, in the present study, smaller moving windows were required, given that small-scale disturbances were analyzed. In addition, it is expected that logging-related disturbances will not be larger in size than 6.5 ha, the minimum mapping unit for the PRODES data set (INPE 2020) used as a forest mask in the present study (see Supplementary Material). Contiguous disturbances larger than 6.5 ha are classified as deforestation, according to the PRODES/ INPE methodology.
Combining composition and configuration in a forest disturbance intensity map A map showing discrete numbers of forest disturbance classes can be derived by combining Pd and Pdd (Zurlini et al. 2006;Riitters et al. 2017). In the following, a justification is given for the use of the chosen metrics to assess disturbance intensity within the study area, before presenting the approach used to combine Pd and Pdd in a forest disturbance intensity map.

Ecological context
Terrestrial disturbances, defined as ''a relatively discrete event causing a change in the physical structure of the environment'' (Clark 1990), vary in many aspects, notably size, spatial distribution, intensity, severity and frequency, which are the components of the disturbance regime (Turner and Gardner 2015). Disturbance intensity is defined as ''the physical energy of the disturbance event per area per time'' and severity as ''the effect of the disturbance event, in organism, community, or ecosystem'' (Turner and Gardner 2015). According to these authors, both components are closely related, since more intense disturbances generally are more severe. An approach to ''measure'' disturbance intensity is to consider the recovery process: the longer the recovery takes, the higher the intensity of the disturbance event (Chazdon 2003).
Disturbances caused by single trees dying and falling (the gap phase dynamics) are the dominant reason for natural forest turnover in tropical forests (Brokaw 1985;Denslow and Hartshorn 1994;Schnitzer et al. 2008;Chazdon 2014). Most canopy gaps created by falling trees are reported to have a gap size ranging from 10 to 100 m 2 (Brokaw 1985;Denslow and Hartshorn 1994;Hunter et al. 2015). Additionally, larger canopy gaps are caused by infrequent Fig. 3 Example of the adjacency (Pdd) metric computation from a binary map. Left: Study area No. 6, binary map of undisturbed/disturbed forests; Center: Moving window (7 9 7 pixels) surrounding the focal pixel (in blue); Right: Adjacency (Pdd) map, containing pixels with the proportion of the surrounding landscape with disturbed-disturbed adjacency. In the case of the focal pixel exemplified Pdd = 77% windstorms or landslides (Nelson et al. 1994;Espírito-Santo et al. 2014;Negrón-Juárez et al. 2017). Artificial gaps created by selective logging activities, despite being designed to emulate natural systems, result in much larger disturbance areas than those caused by the natural processes described above (Chazdon 2014). The amount and the spatial configuration of canopy openings caused by selective logging activities in a given area can be used as a proxy for measuring disturbance intensity.
Composition and configuration were combined in a forest disturbance intensity map comprising five discrete intensity classes: very low (Class I1), low (Class I2), moderate (Class I3), high (Class I4) and very high (Class I5). The reasons for adopting this classification are justified as follows: (i) Isolated pixels in the binary map have a high probability of being just noise (Langner et al. 2018) and even if they correctly represent a disturbance, it can be considered a relatively small one (1 pixel represents an area of 100 m 2 ). The metrics employed here, particularly the contagion metric (Pdd), can be useful for identifying such pixels. (ii) The centers of large disturbed patches are likely to experience different physical conditions than those of small patches (Turner and Gardner 2015). By calculating Pd and Pdd, it is possible to assign different disturbance intensities to the core and edges of the same large forest gap. The core area can be considered more severely disturbed than its edge. Ecologically, this assumption is supported by the fact that a large gap has a very different climate at its center than at the edges, which includes differences in air movement and cooling processes (Oliver and Larson 1996). Resources, such as light availability and soil nutrients, also vary strongly from the edge to the center of gaps (Schnitzer et al. 2008). Furthermore, tree regeneration may be slower at the center of disturbed sites in comparison to at its edges (Chazdon 2014). (iii) The same categorical classes (or types) of forest disturbances, also named ''disturbance profiles'', should have higher inter-disturbance similarity than when comparing different disturbance classes (Lundquist 1995;Zurlini et al. 2006;Buma et al. 2017).
In this study, these categorical classes (k = 5) were constrained by the typical number of logging activities/infrastructure usually found in selectively logged areas (log landings, felling gaps, logging roads and skid trails).

Clustering analysis
Cluster analysis is a common approach used for aggregating information on land cover into meaningful landscape classes. Here, Clustering Large Applications (CLARA), a variant of the k-medoids algorithm developed for the purpose of analyzing large data sets (Kaufman and Rousseeuw 1990), was used to generate the final disturbance intensity map. Clustering was carried out for the aggregated values (the means) of composition and configuration across the range of scales (Fig. 4). The map, aggregated by the mean, is hereafter called a ''multi-scale'' map.
Composition and configuration were combined through CLARA algorithm for k = 5 (five disturbance intensity classes) using the cluster package (Maechler et al. 2018) developed for the R environment (R Core Team 2019).
Assessing the reliability of the forest disturbance mapping Two approaches were adopted to assess whether the categorical classes present in the forest disturbance map correspond to actual forest disturbances. The first one was to compare field data collected in the area with the classes of disturbance. Field data regarding logging activities were collected in 7 of the 10 FMUs under analysis during the dry season of 2017, as a part of the data validation of a previous remote sensing survey (Lima et al. 2019). During the field campaign, GPS location of log landings, felling gaps, logging roads and skid trails were collected. Logging infrastructure (GPS points) that were correctly identified in the maps were included in the comparisons. The second strategy consisted of assessing the rate of forest regrowth amongst different forest disturbance classes. In the section ''Ecological context'', it was stated that an area could have its impact (disturbance intensity) ''measured'' by the rate of its forest regeneration. This ''post-harvest forest recovery'' (White et al. 2019) was assessed using the ''self-referenced normalized burn ratio (rNBR)'' index for each one of the years before logging (2015-2016) and the years after logging (2018 and 2019). A more heavily impacted area will require more time to return to the original values of rNBR than an area that experienced a lower impact. The rNBR is built with bands of the electromagnetic spectrum which are very sensitive to changes in green vegetation and areas with exposed soil (Langner et al. 2016). Pixels eventually disturbed again in the years 2018 and 2019 were excluded from the analysis. For instance, the edges of some FMUs burned in the fire season of the year 2019; therefore, they were not considered in this analysis because they would not be representative for forest regrowth areas.

Composition and configuration
Both Pd and Pdd showed higher values at smaller scales ( Figures S1-S2, Supplementary Material). However, Pdd values stabilized quickly and stayed consistent for landscape extents larger than 5 9 5. In addition, mean values of Pdd did not show great variation among different moving window sizes, ranging from 48.5% at the 3 9 3 and 39.2% at the largest scale (25 9 25). On the other hand, by analyzing Pd values at 3 9 3 extent, the recorded mean was 53.4%, whereas at the largest scale (25 9 25) this value was just 10.4%. The variability in the data, represented by the standard deviation, changed with landscape extent: at the smaller scales, the variance was higher than at larger spatial scales, for both Pd and Pdd. Figure S2 displays the mean values for the different FMUs analyzed. Their patterns are similar: at smaller window sizes, the mean values of Pd and Pdd are high and then decrease as the moving window size increases. Areas No. 6 and No. 8 show the highest Pd and Pdd values, which is an indication that these two areas have experienced (proportionally) larger disturbances in comparison with the other FMUs. The same pattern is observed for the FMUs in the analysis of the final disturbance map derived from the cluster analysis (see more details in the ''Disturbance profiles: analysis of the disturbance intensity classes'' section).  Figure 5 shows the distribution of the individual data points in a scatterplot (represented by each one of the 25,964 pixels from the original land cover map) representing a ''pattern space'' (Riitters et al. 2002). From the smaller window sizes up to the 9 9 9 extension, the metrics have a similar diagonal distribution. However, at 25 9 25 moving window size, most of the individual data points are located at lower values of Pd (\ 40%) and at intermediate values of Pdd, from 30 to 70% (Fig. 5), indicating a decrease in the values of both metrics. Using the mean values of Pd and Pdd (the multi-scale) in the scatterplot creates a ''smoother'' cloud of data points, as the values are no longer categorized following discrete classes. The set of observations at the multi-scale is similar to that obtained when using a 9 9 9 moving window size. The multi-scale values of Pd and Pdd were used to generate a map of disturbance profiles, referred to here as ''disturbance intensity map'', and described in the following section.
Disturbance profiles: analysis of the disturbance intensity classes Mean disturbance proportion (Pd) and mean disturbance adjacency (Pdd) for moving windows with different sizes and for each one of the clusters, derived from the CLARA analysis, are shown in Figure S3. The distinct classes of disturbance trajectories (disturbance profiles) across all landscape extents are shown in Figure S4. Class I1 includes pixels with low values of Pd for all windows sizes. This class covers 23.8% of the disturbance map, considering all FMUs. Figure 6 shows that most of these pixels can be considered as noise, being spatially distant from more clumped spatial configurations. On the other hand, Class I5 includes pixels with the highest values of Pd, ranging from a mean value of almost 90% and decreasing monotonically to a value of 20%, at the larger spatial scale (25 9 25) (see Figure S3). This class represents just 13.9% of all disturbed pixels, and it (mostly) represents the core area of clumped logging-related disturbances, such as log landings and large felling gaps (Fig. 6). Class I4, referred to here as a ''high intensity'' type of disturbance, covers 14.3% of the classified area. It has high values of Pd and Pdd, but they are lower than Class I5 values. In the map (Fig. 6), Class I4 can be associated with the edges of large disturbed areas (log landings and felling gaps) and with logging roads (Fig. 6). The other categories (Class I2 and I3) cover 38% percent of the map and represent ''intermediate'' types of disturbances that can be associated with small felling gaps. However, the distinction between these two classes is not as clear as it is for the other classes (particularly I1 and I5). Figure S5 shows the distribution of disturbance classes for the 10 FMUs. It can be seen that the different FMUs have a diverse logging-related spatial pattern, once the proportional distribution of the disturbance classes does not follow the same pattern for all FMUs. From Fig. 6, this observation is even more evident. Disturbed forests represent 5.1% and 5.5% of the total area of forests of FMUs No. 3 and No. 8, respectively. This value is very similar if the proportion of an area was to be taken as an indicator. However, considering how the proportion of each disturbance class is distributed for both FMUs, it is possible to see that these areas are very different in relation to the spatial pattern. FMU No. 8 has the second highest proportion of Class I5 (very high disturbance intensity), even though it does not have the highest amount of disturbed forests, neither in relation to absolute values nor proportional values (Table S2).
Assessing the reliability of the forest disturbance mapping A total of 165 field points were collected in seven FMUs during the field work carried out in October 2017. These field points represented different types of logging infrastructure: log landings (40 points), felling gaps (61), logging roads (32) and skid trails (32). Not all field points collected within the FMUs could be associated with disturbances on the disturbance map. This is especially true for ''below canopy'' disturbance types, such as skid trails and logging roads. Therefore, from the total of 165 points collected in the field, just 76 points were represented on the forest disturbance map (log landings = 38, felling gaps = 28, logging roads = 8 and skid trails = 2). All other points lie within areas classified as undisturbed forests. Therefore, these 76 points were associated with the forest disturbance classes (Fig. 7). Considering the disturbance classes and the logging infrastructure, log landings can be associated with high and very high disturbance intensity classes, whereas intermediate classes (I2 and I3) were mostly associated with felling gaps. Figure 8 shows the post-harvest forest recovery. The classes defined by the forest disturbance intensity map are clearly distinct in terms of the rNBR values for 2017, the year where the selective logging activities when recorded. This pattern is still recognizable for the year 2018; however, with classes I1 and I2 having almost the same value of rNBR. For the second year after logging, classes I4 and I5, with high and very high disturbance intensities, still show the highest values of rNBR. However, the other classes present virtually the same rNBR values.

Discussion
The main goal of this study was to produce a forest disturbance intensity map to improve the evaluation and monitoring of selective logging activities in areas licensed for timber harvest in the Brazilian Amazon. The main premise of this study is that selective logging activities in tropical forests show heterogeneous patterns and create different types of forest disturbance patches, with different disturbance intensities. Therefore, a well-known methodology designed for large scale studies (Zurlini et al. 2006;Riitters et al. 2017) was adapted to produce a forest disturbance Fig. 7 Field data (logging infrastructure) distribution among different disturbance intensity classes intensity map that accounts for the intrinsic heterogeneity associated with selective logging activities.
The spatial patterns of logging-related disturbances were evaluated throughout multiple observation scales, with the objective of producing and analyzing a map of disturbance profiles (''disturbance intensities''). Two fundamental landscape metrics (Pd and Pdd) were calculated for such a purpose, using different spatial extents. Based on the multi-scale (average) values of Pd and Pdd, a map containing discrete categories of forest disturbance was produced using a cluster analysis (Fig. 6). By analyzing this ''disturbance intensity map'' it was possible to evaluate selective logging beyond counting pixels in a binary (disturbed/undisturbed) map. This is a step further because, in the categorized disturbance intensity map, pixels are classified into different classes of disturbance intensities. In simple binary maps, all logging areas have the same intensity class, based on the mapping of just two land cover classes: disturbed forest and undisturbed forest.

Composition and configuration variability across landscape extents
Composition and configuration, expressed as proportion (Pd) and adjacency (Pdd), respectively, describe spatial patterns of forest disturbance events in numerous landscape scale studies (Zurlini et al. 2006;Zaccarelli et al. 2008;Bourbonnais et al. 2017;Buma et al. 2017;Riitters et al. 2017). These previous studies used spatial pattern analyses to assess disturbance trajectories among different types of forest disturbances (e.g. forest fire, tree harvest, land use change), and clearly distinguish between regional and continental scales. The main goal of the present study was to assess the variation within a specific kind of forest disturbance: tree harvest through selective logging.
Using areas of small-scale and legal (governmentauthorized) logging, it was possible to distinguish between the different types of forest disturbances caused by tree harvest activities, referred to here as forest disturbance profiles. This methodology was originally developed for regional/continental scale studies (Zurlini et al. 2006;Riitters et al. 2017), covering a variety of disturbance types. However, efforts to downscale it have rarely been made (Wickham and Riitters 2019). Here, we have shown that the method can also be used at smaller scales.
The values of Pd and Pdd across different landscape extents showed that selective logging in the study area is composed of disturbances that are locally severe and which are spread throughout a region with few disturbances. This pattern has been described before for forested areas, also dominated by the gap phase dynamics , and for highly anthropogenically-modified landscapes (Zurlini et al. 2007). Similarly, in this study, mean values of Pd are higher and with high variability (represented by the standard deviation) at the smaller scale (3 9 3), whereas at larger landscape extents the standard deviation tends to decrease (Fig. S1).
Using simulated landscapes, Riitters (2019) used the pattern metric space to infer aspects of landscape patterns derived from Pd and Pdd. It is possible to use this theoretical model to interpret the patterns observed in Fig. 5. At small landscape extents, e.g. the 5 9 5 moving window, the cloud of data points in the upper right corner of the pattern space defines focal class pixels characterized as ''clumped'' and defines a landscape with larger and more distinct ''perforations'' or ''holes'', related to the type of background/foreground map under analysis. As we move towards bigger landscape sizes (up to 25 9 25 window size), pixels in the upper right corner of the pattern space are absent. Most of the data points are located at the lower left corner of the metric space (Fig. 5). At this landscape extent most of the observations are classified as ''patchy'', ranging from many small patches to fewer larger patches (Riitters 2019). In practical terms, at this landscape size, most of the observed values of Pd and Pdd are quite similar, regardless of the different assigned disturbance types.
Logging-related disturbances, as observed in the study area, are represented by a mosaic of distinct types of disturbances left after the tree harvest process (Asner et al. 2002;Chazdon 2003;Meijaard et al. 2005;Asner et al. 2009;Chazdon 2014;de Carvalho et al. 2017). These forest disturbances have a different spatial configuration. Recognizing spatial patterns within a given area will depend on the observation scale. The problem of changing observation scales is a long standing issue in landscape ecology (e.g. Levin 1992). This specific author pointed out that, in homogeneous environments, variability (of the system property) will decrease with increasing landscape sizes (moving window sizes). These changes in variability are also seen here across different disturbance profiles: as the moving window size increases, the discrepancy in mean values of Pd and Pdd decreases (Figs. S3 and S4).

Disturbance intensity map: interpretation of disturbance profiles
Unlike previous studies, the approach applied here does not take into account a regional scale, which could contain more than one type of forest disturbance (e.g. forest fires, windstorms, insect outbreaks, clear cutting, shifting cultivation, etc.). Instead it was limited to the boundaries of selected FMUs, covered largely by undisturbed forest (or not directly disturbed forest), given the small-scale nature of tree harvesting activities. The assessment was done at the focal class level (Wu 2004), meaning that only the disturbance class in a binary map was under investigation. By restricting the analysis to FMUs and just focusing on the disturbance class, it was possible to group individual disturbed pixels into meaningful classes of disturbance profiles.
Disturbance profiles have been described as distinct types of disturbances grouped according to their similarity, taking into account certain landscape metrics (Zurlini et al. 2007;Zaccarelli et al. 2008;Buma et al. 2017;Riitters et al. 2017). These profiles may or may not be related to actual, real disturbances types, i.e. a land use or land cover category. For instance, while analyzing the regional patterns of disturbance in the south of Italy, Zurlini et al. (2007) observed that some of the disturbance profiles were correlated with actual land use types, while others were not. This was also observed in the study area.
In the Brazilian Amazon, logging-related disturbances have traditionally been classified in log landings (or logging decks), logging roads (or secondary roads), felling gaps and skid trails (Asner et al. 2002;de Carvalho et al. 2017;Lima et al. 2019). This was the main reason for restricting the cluster analysis to a few classes, which differed from previous studies (Zurlini et al. 2006;Riitters et al. 2017;Wardlaw et al. 2018), and also for not running an algorithm to get an ''optimal'' numbers of k-classes . Considering the FMUs under investigation and the aforementioned logging-related disturbance classes, disturbances caused by skid trails do not result in a large enough canopy cover removal to be captured by optical sensors (Masiliūnas 2017;Lima et al. 2019). However, the input data set captures and successfully maps most of the disturbances that result in direct canopy cover removal (log landings and felling gaps).
Through a visual analysis of Fig. 6, some of the clusters can be associated with actual disturbance classes. Class I5 can be associated with log landings, and Class I4 with logging roads. Class I1 represents mainly noise (thus not real forest disturbances), as it appears in a sporadic, isolated manner, a spatial pattern atypical for selective logging activities. Field points identified as log landings are represented mostly as Classes I4 and I5. These two classes and Classes I2 and I3 can be also be associated with felling gaps (Fig. 7), reflecting the spatial heterogeneity of this kind of disturbance among the FMUs under analysis. Logging roads in this region are difficult to map even with very high resolution data (Lima et al. 2019) while skid trails are basically undetectable with optical remote sensing techniques (Masiliūnas 2017;Lima et al. 2019). In consequence, only disturbance intensity Classes I5 and I4 can be associated with logging infrastructure such as log landings and large felling gaps, on the basis of the collected field points. Some of the field points taken in (probably smaller) felling gaps were mapped as Class I2, which constitutes an intermediate forest disturbance class.
The correlation of forest disturbance profiles with actual disturbances was reported to be more obvious for extreme classes of disturbance intensities (e.g. I1 and I5), considering different disturbance profiles (Zurlini et al. 2007). According to these authors, intermediate classes were more difficult to correlate with real disturbances, which is also observed in our areas. Therefore, while it is possible, at a certain level, to associate the disturbance profiles with logging infrastructure, the forest disturbance intensity map does not perfectly distinguish categories of logging infrastructure. However, with this method it was possible to assign degrees of disturbance intensity, which can have important effects on the regeneration potential of forest stands after selective logging. More heavily impacted areas may experience slower forest succession (Chazdon 2014) and a proliferation of noncommercial pioneer tree species (de Carvalho et al. 2017). This situation is not ideal from an ecological nor from an economic point of view.
Both values of the rNBR, in the year of logging (2017) and in the first year after logging (2018), show a distinct pattern among the forest disturbance intensity classes (Fig. 8). The post-harvest values of rNBR (forest recovery) are higher for more heavily impacted areas, potentially lengthening the time period needed to return to the original, pre-disturbance rNBR values. This is evident for Classes I4 and I5 in year the 2019. The areas of disturbed forest recovered quite fast (Fig. 8), confirming the findings of previous studies carried out in tropical forests (Verhegghen et al. 2015;Dalagnol et al. 2019;Pinage et al. 2019). Furthermore, selectively logged areas classified in the low to moderate disturbance classes (I1, I2 and I3) recovered even more rapidly in comparison with those classified in the high (I4) and very high (I5) disturbance classes. This fast canopy closure after disturbance, detected from optical remote sensing analysis, is in fact caused by the fast regrowth of pioneer plant species that have high photosynthetic activity (Asner et al. 2009;Kleinschroth et al. 2015). In a study carried out in a forest concession in a nearby region (Jamari National Forest, Rondônia, Brazil), Dalagnol et al. (2019) found a similar trend analyzing forest regrowth in function of canopy gap sizes. The authors used aerial LiDAR data, which can capture characteristics of forest structure, showing that large forest canopy gaps will need more time to fully close in comparison with medium or small gaps.
Even though canopy closure after selective logging in tropical forest can occur very fast, the ecological consequences can persist for a long period (Asner et al. 2009;Chazdon 2014). These ecological consequences vary among different types of selective logging infrastructure (de Carvalho et al. 2017), highlighting the importance of mapping categories of disturbance intensities. Therefore, FMUs with high proportions of very high disturbance intensities (Class I5) may not be ideal, taking into account the principles of sustainable forest management. Brazilian legislation presents restrictions on the amount of area, within an SFM, that can be opened up into log landings and logging roads during the harvest process (CONAMA 2009;CEMAAM 2013). However, there are no constraints regarding canopy cover removal for felling merchantable trees, and no further guidelines exist in relation to their spatial distribution. FMU No. 6, for instance has a large proportion of areas classified as I5 due to the spatial distribution of felled trees. In some areas within this FMU, three to four large trees were harvested in close proximity to each other (observed during field work), which caused the clumped spatial pattern evidenced in the disturbance intensity map ( Figure S6). This example shows how forest legislation could additionally consider canopy openings and their overall intensity, as well as their spatial arrangements, to strengthen principles of sustainability in the Brazilian Amazon.

Conclusions
This research had the main objective of mapping forest disturbance intensities in areas under selective logging activities within a focus region in the Brazilian Amazon. These areas were under an SFM regime and had a small-scale type of forest disturbance, which make them difficult to detect and map. We successfully mapped and classified the areas of logging, assigning them classes of forest disturbance intensity and this could help to improve the monitoring of areas under SFM. However, our research has some limitations and its caveats need to be addressed. This study was carried out over 10 FMUs, which in theory should have the same forest management regime. The distinction made between different classes of forest disturbances derived from logging activities was possible because the analysis was constrained to one type of forest related disturbance (i.e. selective logging). The areas of selective logging needed to be clearly separated from other types of disturbances (e.g. forest fire, shifting cultivation, etc.) in order to derive detailed, specific logging-related information on forest disturbance intensities. The approach tested here, despite being successful in separating five distinct forest disturbance classes, can only partially identify logging infrastructure (log landings, felling gaps and logging roads), due to the limitations of the optical remote sensing input data. Further research is needed in this field, e.g. the application of the analysis over areas with a more distinct typology of selective logging infrastructure, as it appears in large forest concessions (Verhegghen et al. 2015). In addition, the post-harvest forest recovery analysis was carried out using exclusively optical remote sensing information. Future research should capitalize on the additional information provided by active remote sensing technologies such as radar and LiDAR.