Remotely-Sensed Carbon Density, Plant Species Diversity and Its Relation to Forest Clearing Dynamics in the Puuc Biocultural State Reserve, Mexico

1. Introduction

It is estimated that about 30% of terrestrial species have been globally threatened or driven to extinction since the 13th century; a process caused by habitat destruction resulting from the expansion of anthropogenic land uses and ecosystem disturbances (Isbell et al., 2022). Novacek and Cleland (2001) estimate that the current rate of massive conversion and degradation of ecosystems by humans further threatens with extinction up to 30% of all extant species by the mid-21st century, an event comparable to some of the mass extinctions that occurred in the past from climatic or geological events. The transformation of ecosystems into anthropogenic land uses is the result of an increasing demand for productive lands all around the globe driven in part by national and international policies for the reduction of poverty and food security. In the face of this ongoing loss, our cultural response has been to increase the number of protected areas to set aside remnants of biological diversity from human dominated areas.

More than 200,000 protected areas exist today, and although the number of protected areas has been increasing globally, these protect only 15% of the planet’s lands and between 10 and 14% of the planet’s water (UNEP-WCMC and IUCN, 2016). Tropical forest reserves are especially important in this context since they preserve 68% of the global carbon stored in terrestrial ecosystems and harbor the highest degree of biological diversity in the planet (Bebber and Butt, 2017; Laurance et al. 2012).

Overall, protected areas have been effective. Bruner et al. (2001) suggest that the majority of parks in tropical countries have been successful at stopping land clearing and mitigating logging, hunting, fire, and grazing. A national level study by Blackman et al. (2015) in Mexico, as well as global studies by Laurance et al. (2012) and Coetzee et al. (2014) provide evidence that parks have preserved species and their tropical forest habitats.

Nonetheless, tropical protected areas are often criticized for not providing adequate protection to biodiversity. Data published by Laurance et al. (2012) reveals that about half of all reserves have been effective or performed passably, but the rest are experiencing an erosion of biodiversity that is often alarmingly widespread taxonomically and functionally. As stated by Rodriguez and Rodriguez-Clark (2001), underfunded and understaffed park agencies are challenged by pressure exerted by human populations living around and within national parks and other protected areas. The fact is that these protected areas, once isolated from highly populated areas and distant from threats, are now embedded in social-ecological land systems characterized by patches of preserved or managed natural ecosystems in a ‘matrix’ of urban, agriculture and livestock ranching that has expanded rapidly, increasing land conflicts between stakeholders as well as degradation outside and inside the protected areas.

In the case of the Puuc Biocultural State Reserve (PBSR), a tropical dry forest dominated reserve located in the Central Yucatan Peninsula, Mexico, effectiveness depends on the coordination of multiple communities and entities that form this special state park. The PBSR was established in 2011 by the Yucatan State Government primarily as a means to preserve the biological diversity of the region, its historical value for Mayan culture and the sustainable livelihoods persisting in the region (Secretaria de Desarrollo Sustentable, 2022). This type of reserve is the first of its kind in Mexico. It covers the territory of five municipalities that have historically been inhabited by low density human populations dedicated to small traditional agricultural practices (milpas) and cattle ranching in the form of communal properties (ejidos) and private properties. Some of these communities implement sustainable forestry projects as well. The PBSR was established as a corridor connecting the forested areas of the south of Campeche and Quintana Roo states with the Celestún Biosphere Reserve in north-eastern Yucatan.

The main activities that cause forest disturbances in the PBSR are agricultural activities such as cattle ranching (which is concentrated in small areas in the north and southeast of the reserve), the cultivation of vegetables such as corn, beans, sorghum and citric fruits, and forestry activities related to charcoal production and house construction. These activities are mostly implemented in milpa systems. Milpa systems are traditional Mayan land management schemes that involve slash-and-burn agriculture of multiple crops. Crops are followed by fallow after at least 4 years which allows for forest regrowth and soil nutrient recovery (a period called ‘barbecho’).

The PBSR management plan requires community-based planning of all productive activities coordinated by an inter-municipal board. It is important to note that 55% of the reserve is managed by ejidos while the rest are private properties. All activities are planned to minimize environmental impact and safeguard the ecological sustainability needed to preserve the ecosystem services that the region provides to its inhabitants and the surrounding communities (Secretaria de Desarrollo Sustentable, 2022). Before its creation, more than 86,000 hectares had been reported as forested areas within the reserve which represents more than 63% of its extent. The PBSR strives to maintain a tight balance between productive activities and forest conservation through community-based governance. But the intermunicipal board has also identified threats to this balance. The PBSR has documented forest loss and disturbance through the expansion of cattle ranching (often incentivized by government subsidies); unsustainable forestry and hunting activities; road construction and increased forest clearing due to the shortening of the barbecho or fallow period. The extent to which these activities have transformed the landscape after the creation of the PBSR has not been assessed. More than 10 years after its creation, an assessment of the forest cover dynamics before and after the creation of the reserve can help understand the role of the PBSR in this socioecological system.

Over the past several years there have been significant advances in the design of continuous land cover change (CLCC) mapping algorithms that take advantage of the high-quality Landsat data archive that became freely available in 2008 (Cohen et al. 2017; Zhu, 2017). CLCC algorithms produce annual land cover change maps showing the location and extent of deforestation events occurred in every single year. In addition, numerous applications have used the power of multispectral and hyperspectral optical, radar, and LiDAR sensors for mapping ecosystem age, biomass and carbon stocks, alpha and beta diversity, biochemical heterogeneity, tree height, vertical structure, and complexity, among other attributes (Wang and Gamon, 2019). In tropical forests, state-of-the-art methodologies using multispectral, hyperspectral imagery and/or LiDAR usually achieve R2 values of ~0.5–0.7 for canopy functional trait mapping, while categorical maps of plant functional types can achieve overall accuracies of ~0.8 (Portillo-Quintero et al. 2021).

In this study, we integrate state-of-the-art remote sensing tools and products for assessing forest clearing dynamics inside and around the PBSR. We combined spatially-explicit models of carbon density and tree species richness, derived by using L-band radar backscatter data as well as radar texture metrics, climatic and field data (Hernandez-Stefanoni et al. 2021); and an original forest loss product derived from the analysis of time series of Landsat 5, 7 and 8 from the year 2000 until 2021. In the context of the PBSR management objectives, the goal of the study is to answer the following questions: a) how has the creation of the PBSR impacted the frequency and area of forest loss within its boundaries? b) Are productive activities contemplated in the PBSR governance model compatible with the maintenance of carbon density and species richness?.

2. Methods

2.1. Study Area

The study area was defined as the area within the PBSR boundaries and the adjacent and wider landscape around the reserve (Figure 1). This area was defined by all municipalities intercepting a 25 km area of influence around the PBSR. This included municipalities immediately adjacent to the PBSR and municipalities adjacent to these. Some municipalities extend far beyond the 25 km area of influence and therefore were only represented partially in the analysis. Municipalities completely included in the study area were Oxkutzcab, Santa Elena, Muna, Ticul, Akil and Dzan in the Yucatan state. The study area included partial coverage of the municipalities Hopelchen, Calkini and Hecelchakan in Campeche state and the municipalities of Tekax, Tzucacab, Chacsinkin, Tixmehuac, Teabo, Mani, Mama, Chapab, Sacalum, Abala, Opichen, Kopoma, Maxcanu and Halacho in the Yucatan state.

The vegetation communities in this region are predominantly secondary semi-deciduous forests that can reach a canopy tree height of 15-20 m, losing 50%-75% of its leaves during the dry season. Annual rainfall ranges from 1000 to 1500 mm with an annual temperature of 25.9C to 26.6C. Elevation of the terrain in the area ranges from 60 to 250 m above sea level. The dominant species are Neomillspaughia emarginata, Gymnopodium floribundum, Bursera simaruba, Piscidia piscipula and Lysiloma latisiliquum. The local geomorphology combines flat areas with limestone hills of moderate slopes. The landscape consists mainly of forest land (94% of the area), mostly dominated by forest stands of different successional ages.

2.2. Semi-automated forest disturbance detection from 2000 to 2021.

2.2.1. Generation of a baseline map

The detection and mapping of forest loss events was performed using a semi-automated continuous change detection method. A semi-automated method was necessary in order to avoid overestimation of loss due to the inclusion of errors related to changes in biomass in agricultural crops (Portillo-Quintero & Smith, 2018). Overall, the method consisted in detecting changes in vegetation index values every three months over known forested areas (post-classification change detection) from the year 2000 until the end of the year 2021.

First, we produced a baseline forest cover map for an initial time step in the year 2000. We downloaded a cloud free (<5%) Landsat 5 (L5) scene (20/46) surface reflectance product for a single date (March 28th) for the year 2000 using the USGS EROS Science Processing Architecture (ESPA) online interface. Then, we used R packages (raster, rpart, randomForest, gdal and rgeos) for image pre-processing, cloud masking and supervised classification using the Random Forest algorithm.

For classification, we collected training data through on-screen interpretation of the L5 scene in a RGB453 false color composite and very high-resolution imagery available in ArcGIS Online (ESRI), in combination with the INEGI VI Series Land Cover Classification map. We arbitrarily chose a quantity of ~800 locations as the target for training data collection. The number of training polygons per class was calculated through proportional allocation using the original land cover proportions in the INEGI land cover map. The proportional allocation equation used was Sh=(Nh/N)*n where Sh is the sample size per class, Nh is the class size, N is the total landscape size and n is the desired total number of samples. A total of 836 polygons were digitized with 595 polygons for the forest class and 241 polygons for the Non-Forest class. The forest class represented all forest types defined in the INEGI land cover map, and the Non-Forest Class represented anthropogenic land covers and water features. Finally, the Random Forest classifier was performed using visible, near infrared and short-wave infrared bands as well as a Normalized Difference Moisture Index (NDMI) and a Normalized Burn Ratio 2 (NBR2) products from the L5 scene.

To evaluate the accuracy of the baseline map, we implemented a stratified random sampling approach for generating validation points and an area-based error matrix (Congalton, 1991; Oloffson et al. 2014). Through this method we calculated overall, user’s and producer’s accuracy with error-adjusted area estimates calculated at a 95% confidence interval. Based on the extent of Forest/Non-forest land covers, we estimated the optimal size of samples per class for accuracy assessment. We used the proportional allocation formula to estimate the number of samples per class. We targeted a total of 800 samples for validation purposes as a match to the number of samples used for training the image classification. Proportional allocation formula resulted in 621 samples for the forest cover and 179 samples for the non-forest cover. Samples were collected using a 250 x 250 m grid over the baseline map in ArcGIS 10.7. Cell grids per class were randomly selected and a centroid point was generated for each cell. These centroid points were used for map-to-reference data comparisons using Google Earth Pro.

2.2.2. Post-classification change detection

For post-classification change detection, we used the Normalized Burn Ratio 2 (NBR2) vegetation index products from Landsat 5, Landsat 7 and Landsat 8 SR tier 1 image collections. We used the NBR2 given its overall better performance than other indices in change detection algorithms (Smith et al. 2019). Using Google Earth Engine (GEE), we filtered Landsat imagery to four trimesters in every year from January to March, 31st (Q1), April 1st to June 30th(Q2), July 7th to September 30th (Q3) and October 1st to December 31st(Q4). For each trimester, we selected only clear observations, masking clouds and cloud shadows to every scene using the pixel qa product; and generated a composite of maximum NBR2 values for every trimester.

Each NBR2 trimester composite was inspected in ArcGIS, and a threshold for values that were <-1.5 standard deviations from the mean was applied in order to isolate non- or little vegetated surfaces. This process was done interactively through visual inspection and manual thresholding in ArcGIS. In some cases, thresholds were slightly adjusted to values between <-0.5 and <-1.5 standard deviations to match the exact dimensions of non-vegetated surfaces. Binary maps of vegetated (0) and non-vegetated surfaces (1) were generated.

We subtracted these rasters from the forest mask of the previous year in order to isolate forest disturbances. We then added all trimester forest disturbance masks in order to produce annual forest disturbance maps. Following Hamunyela et al. (2020), we flagged forest disturbances only after two or more consecutive observations of low NBR2 values, however, in some cases we allowed all possible observations to be included in the annual product. This was done given that cloud contamination and L7 missing data affect the probability of detecting consecutive observations. Changes in trimester NBR2 values following cultivation of recently cleared land can also mask the detection of these consecutive observations. Isolated pixels were successfully filtered using a focal 3x3 majority filter leaving the more regularly shaped forest disturbances intact. For every time step, we achieved a final annual product only after satisfactory visual inspections. Final deforestation polygons where limited to areas larger than 0.18 ha in order to eliminate spurious isolated pixels.

For assessing accuracy of the annual forest clearing product, we implemented a two-stage cluster sampling design to generate validation locations. Using a 5 km x 5 km grid we calculated the percentage of forest clearings detected per cell using the combined annual forest clearing dataset. This allowed us to generate a grid with forest clearing percentages. We categorized the grid using the Jenks optimization method (Natural breaks) in ArcGIS in four classes of stable forest or no loss (128 cells), low forest loss (1,005 cells), intermediate forest loss (194 cells) and high forest loss (29 cells). We then used the Neyman optimal allocation formula (Stehman et al. 2011), which takes into account the number of cells and the standard deviation of forest cleared area per cell, to calculate the number of sampling points per class needed for validation in 25% of the landscape. We targeted 800 samples for validation purposes as a match to the number of samples used for training the image classification and for validation of the baseline map. The Neyman optical allocation method estimated that proper representativeness of the landscape would require selecting 9 cells in the stable forest class, 627 in the low clearing class, 124 in the intermediate clearing class and 40 in the high clearing class. Within these set of cells, we re-calculated the area with stable land use and land cover (no forest cover change) and with forest change (forest clearings) and further estimated the amount of sampling points needed for validation by proportional allocation. This resulted in 593 sampling points generated to cover stable land cover and 207 sampling points to cover the forest-loss class (for a total of 800 validation points). Each of the samples were individually inspected in Google Earth Pro and using the historical imagery tool, we looked for evidence of change between 2000-2021 to validate the map class. An area-weighted confusion matrix was used to compute overall, producer’s and user’s accuracies, along with standard errors and forest clearing areas estimated at a 95% confidence interval (Stehman, 1997; Potapov et al., 2014).

2.3. Carbon density and tree species richness models

The methodology to derive carbon density and species richness maps for the study area is thoroughly explained in Hernandez-Stefanoni et al. (2021). In their study, authors generated spatial models of tree species richness and carbon density for the entire Yucatan Peninsula of Mexico (Figure 2). Field data on tree species richness and aboveground biomass (Mg / ha; used as a measure of carbon density) was obtained from the National Forest Inventory (NFI) of Mexico from surveys conducted between 2009 and 2014. A total of 2847 plots were used and these were distributed in the forested areas in a stratified sampling scheme depending on the vegetation type. Each field plot consisted of 4 circular 400 m2 sub-plots distributed within an area of 1 ha. All trees with a diameter at breast height (DBH, 1.3 m) ≥ 7.5 cm were identified at species level and their diameters and heights were measured. AGB stocks per plot were calculated in standard units (Mg ha−1) using local and regional allometric equations, which were developed by vegetation type, plant growth form and DBH size.

The study used ALOS-2 PALSAR-2 (Advanced Land Observing Satellite Phased Array L-band Synthetic Aperture Radar) mosaic tiles for the year 2015—with 25 m resolution and dual polarization (HH and HV) preprocessed for orthorectification, slope correction, and radiometric calibration. It used backscatter coefficients (γ°) and Normalized Difference Backscatter Index (NDBI) products as well as variability in backscatter values among pixels using texture metrics and climatic water deficit (CWD) data. This process generated 24 variables, considering 8 texture metrics and three bands (HH, HV and NDBI) that were used for AGB and species richness estimations. A random forest regression model was used to estimate the spatial distribution of AGB and tree species richness from all explanatory variables. Model validation of the regional maps obtained herein using an independent data set of plots resulted in a coefficient of determination (R2) of 0.28 and 0.31 and a relative mean square error of 38.5% and 33.0% for aboveground biomass and species richness, respectively, at pixel level.

2.4. Statistical analysis

In order to assess whether the creation of the PBSR has impacted the frequency and area of forest loss within its boundaries, we calculated basic statistics on the number of forest clearing events per year (frequency) and forest loss per year (in hectares/year) in three different geographic analysis units before and after the year 2011 (date of its creation): inside the PBSR boundaries, in a 12.5 km buffer ring around the PBSR boundaries and in a 25 km buffer ring around the PBSR. The two buffer distances were arbitrarily selected in order to compare forest clearing dynamics in the communities within the PBSR and the communities around it. We then performed pairwise comparisons between geographic analysis units before and after 2011 using the Kruskal-Wallis rank sum test (kruskal.test()) in R-software.

Using the Mann-Kendall trend statistic package (Kendall) in R, we also investigated if there were significant trends in the time series of annual forest area and evaluated whether these trends were different before and after the creation of the PBSR.

In order to further investigate if there were significant spatio-temporal patterns in the frequency of forest clearing events before and after the creation of the PBSR, we used the ‘Emergent Hotspot Analysis’ or EHSA tool in ArcGIS Pro (ESRI). The EHSA tool aggregates points into space-time bins. Aggregation of points per year was performed at a 2.5 × 2.5 km bin size in one year intervals using hexagon grid cells. To get a measure of the intensity of feature clustering in each bin, the hotspot analysis uses a space-time implementation of the Getis-Ord Gi* statistic, which considers the value for each bin within the context of the values for neighboring bins in space and in time (ESRI, 2023). We selected only adjacent bins (6) to each hexagon cell as neighboring bins and a temporal frame of 3 time steps (3 years). The time series of Getis-Ord Gi* z-scores for each bin and its spatio-temporal neighbors is then assessed using the Mann-Kendall statistic. The result from this analysis is a clustering trend z-score, p-value, and binning category for each location. Binning categories are based on the intensity of the clustering (Z-scores and p-values) for high number of forest clearing events (hotspots) and low numbers of forest clearing events values (cold spots). Hotspots and coldspots are automatically classified into new, consecutive, intensifying, persistent, diminishing, sporadic and oscillating bin categories based on clustering patterns in time and space.

We then assessed whether the PBSR governance model minimizes the impact of productive activities on forest carbon density and tree species richness. For this, we selected only the annual forest loss detected after 2015 (the year of our reference biomass and diversity map) and calculated the mean carbon density and the mean tree species richness for each forest clearing polygon from 2015 until the year 2021. Using a Kruskal-Wallis rank sum test implemented in R, we tested for differences between the distribution of carbon density and species richness values in forest loss polygons and an equivalent random sample in stable forest areas for two geographic analysis units: within the PBSR and outside of the PBSR up to a surrounding area of 25 km distance from its boundaries.

3. Results

3.1. Forest cover baseline and annual forest clearing dataset

The semi-automated classification of the Landsat 5 scene surface reflectance product for March 28th, 2000 allowed us to produce a forest cover map that was used as a baseline for detecting changes in forest cover at an annual basis (Figure 3). The baseline forest map was evaluated using a confusion matrix which showed a high overall accuracy of 94% with producer’s accuracy of 62.7% for the forest class, and 91.5% for the non-forest class, and a user’s accuracy of 96.8% for the forest class and 82.6% for the non-forest class (Table 1). Because the land cover matrix is diverse and it includes tree plantations and horticultural crops, pixels representing these land uses would likely be included in the forest class and vice versa, forest class pixels representing early successional stages would likely be included in the non-forest class. Nonetheless, we considered the accuracy of the map as satisfactory for delineating the overall forest cover distribution by the year 2000.

The map shows that the total error-adjusted estimated area of forest cover for the year 2000 was 775,846 ha (±15,001 ha at a 95% confidence interval) representing 79.9% of the total terrestrial landmass of the study area which includes the PBSR and a surrounding buffer of 25 km. Forests are distributed evenly across the study area with the largest continuous fragments found near and within the PBSR. Non-forest areas are represented mostly by rural anthropogenic settlements, towns and small agricultural developments (mostly milpa) along highway 184 in the northeastern part of the study area in the state of Yucatan. The northwestern part of the study area, in the state of Campeche, also has large areas of non-forested land represented mostly by intensive agriculture developments and towns along highway 180.

The annual deforestation dataset was successfully generated with an overall accuracy of 91% (Table 1). The matrix shows that the forest clearing class had a moderate rate of commission errors (user’s accuracy of 76.9%) which means areas of no change were counted as forest cover change on a 23.1% rate, while it shows a lower rate of errors of omission (producer’s accuracy of 83%). We considered the results satisfactory for showing overall forest clearing trends spatially and temporally, and comparable to other studies at a similar geographic scale (Diaz-Gallegos et al. 2010, Griffiths et al. 2018).

The estimated error-adjusted area of detected annual forest clearings from the year 2000 until the year 2021 was 230,511 ha in total (±19,979 ha). Results also show a mean and median forest clearing patch area of 1.26 ha and 0.54 ha respectively (Std Dev = 3.36), ranging from patches of 0.20 ha to large forest clearings of 229 ha. The majority of patches (98%) fall within the 0.20-7.27 ha range, with the rest falling mostly under 50 ha and only five patches greater than 100 ha.

3.2. Comparisons by geographical unit and time period

Over the entire study area, annual forest clearings occurred on an average of 1,552 clearings per year with an average area of 1,951 ha per year. The extent of deforestation peaked for the year 2009 with 6,237 ha of forest cleared in 4,272 events. The year with the lowest forest clearing detected was 2001 with only 261 ha of forest cleared. In general, the dynamics of forest removal was highly variable during the whole study period. We observed the highest and the lowest forest area cleared during the first ten years. After 2011, clearings peaked on 2019 with 3,154 ha cleared in 1,948 events. Spatially, clearings were mostly clustered outside the boundaries of the PBSR and less frequent inside the PBSR (Figure 3). Most of the clearing patches greater or equal to 50 ha (n=22) were located in the southeastern and eastern part of the study area inside or near to the state of Campeche and on the northern part of the study area in the state of Yucatan. In Campeche, all of these corresponded to clearings in landscapes dominated by commercial mechanized agriculture while in Yucatan state, these corresponded partially to large contiguous patches of milpa agriculture and some dedicated to commercial mechanized agriculture.

When analyzed by geographical unit (inside the PBSR and in buffer rings of 12.5 km and 25 km around the reserve), forest clearing dynamics greatly differed (Figure 4). Dynamics were very low and stable inside the PBSR during the whole study period (Figure 4A). The average forest clearing area inside the PBSR was 73 ha per year from an average of 60 events per year (for an average 1.21 ha per clearing patch). The maximum area cleared was recorded for 2009 with 280 ha in 180 clearing occurrences. In the 12.5 km buffer, annual forest clearing was higher with an average of 698 ha cleared per year and an average of 582 clearing occurrences per year (for an average 1.19 ha per clearing patch). Similarly, peak forest area cleared was recorded for the year 2009 in the 12,5 km buffer with 2,797 ha cleared on that year. In the 25 km buffer, the magnitude of area cleared and number of occurrences was higher for the years 2000, 2002, 2008, 2009 and 2019 in comparison to the other geographical units of analysis (Figure 4C), with the largest area cleared registered for the year 2000 and followed by the year 2009. Kruskall-Wallis tests shown that forest clearing area and frequency per year were significantly different between geographical units (Table 2). Figure 5 shows the differences in the ranges of forest area and frequency between geographical units before and after the creation of the PBSR in 2011. In general, we observed a general shift in the magnitude of forest cleared before and after the creation of the PBSR, especially outside of the PBSR but the differences were not significant in these two time periods for any of the geographical units (Table 2).

Mann-Kendall trend tests showed no temporal trends detected in the area cleared during the whole study period for the three geographical units (Table 2). There was also no evidence of significant positive or negative temporal trends in forest area clearing dynamics before and after the creation of the PBSR, either inside of the reserve or in the 12,5 km and the 25 km buffer rings (Table 2).

3.3. Emergent hotspot analysis

Emerging trends were assessed using the space-time hotspot analysis in order to investigate whether the creation of the PBSR impacted the spatio-temporal trends in forest clearing. The analysis is presented for forest clearings detected between 2000 and 2011, and between the 2012-2021 period. Figure 6 shows the results of this analysis and depicts the distribution of areas of stable low clearing intensity (e.g., persistent cold spots: low frequencies and clustering maintained in time and space during the period of study) and stable high intensity (e.g., persistent hotspots: high frequency and clustering maintained in time and space during the period of study) of forest loss events. It also shows areas where recent increasing trends and decreasing trends in forest clearing clustering are occurring (e.g., new hotspots and cold spots) and where high clustering happens occasionally (e.g., oscillating hotspots and cold spots). Areas that were not categorized did not show any significant trend across the years.

Results indicate a widespread occurrence of oscillating high intensity areas before 2011 with one cluster of persistent hotspots in the northern part of the study area. Before 2011, newer or recent hotspots were spread out across the study area. Persistent and oscillating low intensity cold spots were evident in the central part of the study area corresponding with the eastern and western section of the PBSR. After the creation of the PBSR in 2011, stable high intensity areas were reduced to smaller ranges around towns in the northern part of the study area and some locations in the southern area that correspond to the state of Campeche. We observe, after 2011, a reduction of the hotspots in the middle of the PBSR and more empty cells all across the PBSR, which means there was no sufficient data to compute the spatio-temporal statistical analysis. There was also a reduction of the number of persistent cold spots inside the PBSR. New hotspots after 2011 are mostly located in the northern part of the study area.

3.4. Comparison to carbon density and species diversity models

Results show that slash-and-burn activities select areas of lower carbon density and species richness in comparison to what is available. This pattern was detected inside and outside the PBSR. Figure 7 shows histogram comparisons of the frequency of forest clearing events by carbon density and species richness ranges. Frequencies are compared between carbon density and species richness in forest clearing polygons after 2015, and a similar number of random locations in forested areas. Inside the PBSR, forest clearing events cluster around the 99-150 Mg/ha carbon density range and the 33-40 species/ha range. Outside the PBSR, forest clearing is more widespread and occurs around 42-152 Mg/ha. It also occurs in a bimodal pattern in relation to tree species richness with a cluster in the ranges of 11-19 species/ha and 27-41 species/ha. Kruskall-Wallis tests show that the distribution of forest clearing events is significantly different from random in all pairwise comparisons (Table 2) and is shifted towards lower values in clearings, especially outside the PBSR (Figure 7).

4. Discussion

The implementation of a semi-automated method for mapping annual forest clearings for the 2000-2021 period using the Landsat archive allowed us to accurately quantify forest cover change dynamics inside and outside the PBSR. However, results from the accuracy assessment suggest that omission errors were high for the forest class in the 2000 baseline map (Producer’s accuracy 62.7%) and for the annual forest clearing dataset (producer’s accuracy 53.9%). This omission error rate might indicate an underestimation of forest clearing activities and suggest that the conversion of forests to milpa and other anthropogenic land uses might be occurring at a higher frequency. The use of a post-classification change detection method using NBR2 thresholding might have underestimated the occurrence of forest clearing in some instances. This could happen if, after moderate to high rainfall events or land irrigation, a trimester quantification of the NBR2 threshold does not accurately reflect the separation between forest cover and anthropogenic land uses. This masking effect could occur given that NBR2 uses short-wave infrared bands highly sensitive to surface water content.

In spite of these potentially omitted clearing events, the dataset registered 34,138 forest clearing events in a pattern consistent with the spatial distribution of slash-and-burn activities expected for landscapes undergoing socioecological processes where a conservation area has been in place. The unbiased (error-adjusted) estimation of total forest area cleared during the 2000-2021 period was 230,511 ha (±19,979 ha) which adjusts estimates based on omission errors and commission errors calculated in the accuracy assessment.

Regarding the impact of the creation of the PBSR on the area and frequency of forest clearings, our results suggest that although forest clearing was significantly more intensive outside of the PBSR than within the PBSR during the entire 2000-2021 period, there is no evidence suggesting that the frequency and magnitude of forest clearing has changed over the years after the creation of the PBSR in 2011. Figure 4 and Figure 5 show slight differences in the magnitude of forest clearing events from 2000 to 2021, but this pattern was evident only for the 25 km buffer ring outside the PBSR. In fact, Mann-Kendall trend statistics suggest there is no trend in the annual forest clearing time series in any of the geographical units (entire study area, inside PBSR, in the 12.5 km buffer ring and in the 25 km buffer ring) and for any of the time periods (pre-2011 and post-2011). The EHSA analysis does show, however, that forest clearing clustering (hotspots) during the 2012-2021 period was less widespread and was confined to areas outside the PBSR, and that forest clearing clustering is on a downward trend within the PBSR (Figure 6). In fact, the EHSA analysis was unable to find sufficient forest clearing data in numerous cells within the boundaries of the PBSR that would yield any significant spatiotemporal trend during the 2012-2021 period. EHSA results therefore suggest a general shift towards the reduction of forest clearing intensity in the study area that could be a product of landscape scale processes not linked to the creation of the PBSR such as migration and socioeconomic processes (Secretaria de Desarrollo Sustentable, 2022) but it also suggests a minimizing effect of the PBSR on the frequency of forest clearing events inside the reserve.

For the 2000-2021 period, published literature and national reports have mixed results regarding the trends in forest clearing for the Yucatan peninsula. Krylov et al. (2018) reports a higher annual rate of deforestation for the Yucatan Peninsula in comparison to national rates. Previous research by Portillo-Quintero & Smith (2018) also lists the Yucatan Peninsula as one of the largest hotspots of forest clearing of the Mesoamerican region. Lawrence et al. (2020) estimated forest losses up to 85% in multiple individually managed and parcelized ejidos in the Yucatan state between 1985 and 2016. Nonetheless, Ellis et al. (2017) indicates that even though annual deforestation is high in the State of Yucatan, it remains lower than in Quintana Roo and Campeche. When observed at the scale of the PBSR and its surroundings, a rapid report using the public online portal of the Global Forest Watch shows that tree cover loss in the same study area (PBSR and its surroundings) does not follow a negative or positive trend and has been stable since the year 2001 until the year 2021. Our results show that forest clearing activities have been reduced and/or remained low and stable within the PBSR and also, forest clearing has remained stable (although less intense according to EHSA) in the surroundings of the PBSR. These results indicate that the study area has experienced lower deforestation rates compared to the reported statewide and regional deforestation processes in the Yucatan peninsula.

The comparison of the annual forest clearing dataset to remotely-sensed carbon density and tree species models further reveals that the PBSR is controlling agriculture dynamics within its boundaries. Our results show that land owners are clearing forests with lower carbon density and tree species diversity in relation to what is available for them at the landscape scale. This is expected given that milpa systems target early to intermediate secondary forests after several years of abandonment (barbecho). However, our data indicates that carbon density and species richness in cleared areas were generally lower (more skewed to the left) relative to random sites outside the PBSR than inside (Figure 7), suggesting that land owners within the PBSR might be practicing longer barbecho periods thereby allowing forests to gain higher carbon density and tree species diversity before slash-and-burn compared to landowners outside the PBSR. Allowing a longer barbecho period is an important component of the traditional shifting cultivation strategy that allows the recovery of soil fertility and results in crops of higher yields (Uribe-Valle & Petit-Aldana, 2007). Bray et al. (2004) reported similar patterns in milpa systems of Central Quintana Roo where longer barbecho periods and the maintenance of late secondary forests are attributed as causes for lower observed deforestation. The PBSR management plan clearly indicates that the practice of short barbecho periods threatens the maintenance of biodiversity and sustainability of agricultural systems in the region and therefore, promotes the use of longer barbecho periods.

5. Conclusion

Since its creation in 2011, the management goal of the PBSR is to preserve sustainable livelihoods for the ‘ejidatarios’ and private land owners that are part of it. It also aims to secure the provision of ecosystems services for a human population of more than 130,000 living in its surroundings. However, dynamics of forest clearing outside of the PBSR have been affected by agricultural intensification (Secretaria de Desarrollo Sustentable, 2022, Proust et al. 2015). The intensification of agriculture threatens to negatively impact the sustainability of the milpa system by shortening the barbecho periods in an effort to maximize productivity. The PBSR is also currently threatened by the expansion of mechanized agriculture that is predominant in the Campeche section of the study area. By combining products from the analysis and transformation of remotely-sensed data (Landsat archive; ALOS-2, PALSAR-2), our study was able to reveal lower forest clearing rates suggesting different underlying socioecological processes occurring in the study area compared to patterns of forest clearing reported at the state and regional scale. Our results suggest that the PBSR effectively acts as a stabilizing forest governance scheme that reduces the impact of productive activities by lowering the frequency of forest clearing events and preserving late secondary forests.

While the integrated analysis of remotely-sensed forest attributes such as biomass, carbon density and species richness requires extensive efforts in data collection at significant financial costs, the implementation of continuous change detection methods using the Landsat or Sentinel archive can be accurately achieved with open-source GIS software and cloud-based platforms at low costs. In this study, we combined carbon density and tree species diversity models at 25 m pixel resolution that were generated by leveraging national scale forest inventory data, with a Landsat-based continuous change detection product. Given the urgency for more integrated assessments that can help guide conservation efforts, the remote sensing community should continue to concentrate efforts in providing alternative optimal field data collection strategies and modeling techniques to spatially predict key tropical forest attributes. These efforts should target the minimization of financial costs while still yield accurate results. Ultimately, the integration of remote sensing based models at sub-national and local scales can speed up the generation of knowledge about the spatial variability in socioecological processes and generate information better adapted to forest governance scales.