Predictive Modeling of Future Forest Cover Change Patterns in Southern Belize

Tropical forests and the biodiversity they contain are declining at an alarming rate throughout the world. Although southern Belize is generally recognized as a highly forested landscape, it is becoming increasingly threatened by unsustainable agricultural practices. Deforestation data allow forest managers to efficiently allocate resources and inform decisions for proper conservation and management. This study utilized satellite imagery to analyze recent forest cover and deforestation in southern Belize to model vulnerability and identify the areas that are the most susceptible to future forest loss. A forest cover change analysis was conducted in Google Earth Engine using a supervised classification of Landsat 8 imagery with ground-truthed land cover points as training data. A multi-layer perceptron neural network model was performed to predict the potential spatial patterns and magnitude of forest loss based on the regional drivers of deforestation. The assessment indicates that the agricultural frontier will continue to expand into recently untouched forests, predicting a decrease from 75.0% mature forest cover in 2016 to 71.9% in 2026. This study represents the most up-to-date assessment of forest cover and the first vulnerability and prediction assessment in southern Belize with immediate applications in conservation planning, monitoring, and management.


Introduction
Throughout the world, deforestation, degradation, and fragmentation threaten the integrity of tropical forests and the biodiversity that they contain.About half of the world's tropical forests have been cleared [1], and between 1980 and 2000, over 80% of new agricultural land originated from forests [2].This deforestation has unknown long-term effects on biodiversity, rare species, ecosystem processes and functions, climate patterns, and the existence of important resources such as medicines and crop relatives [3,4].Considering tropical forests' critical function in sustaining biodiversity, ecosystem services, and local livelihoods, understanding the patterns of forest cover loss and implementing conservation actions in strategic locations to prevent deforestation is crucial.
Belize is generally recognized as a highly forested country with a total mature forest cover of 57.6%, according to unpublished data from February 2017 [5].The country does not necessarily contain "primary" or "virgin" broadleaf forest due to its history of land use practices performed by the ancient Maya and natural and anthropogenic disturbances, such as hurricanes, fires, and logging [6].Therefore, areas that have not been cleared since 1980 are referred to in this study as "mature forests" instead of "primary forests".The Maya Golden Landscape (MGL), located in the Toledo District in southern Belize, is still mostly forested and has retained a greater amount of mature forest cover than other areas in Belize.The 311,610-hectare MGL is a mosaic of private and national protected areas, private lands, and Maya and Hispanic communities (Figure 1).The region includes the primary biological corridor in southern Belize, which is the only remaining broadleaf forest link between the Maya Mountains and the lowland broadleaf forests that extend to the coast.This connection is critically important on both a national and regional scale as part of the Mesoamerican biological corridor.The MGL is also part of Mesoamerica's Selva Maya, which is the second largest remaining tropical rainforest in the Americas, after the Amazon.The Selva Maya supports over 400 bird species, about one billion overwintering migratory bird individuals, and critical populations of threatened species such as the jaguar, Yucatán black howler monkey, and Baird's tapir [6].
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 17 typically fine scale and site-specific [33].This method could be applied more widely to conservation planning in other landscapes, especially regions with severe financial constraints, as a way to prioritize the protection of threatened areas.

Forest Cover Classification and Change Analysis
Forest cover was classified utilizing a supervised classification of Landsat 8 Operational Land Imager (OLI) data in Google Earth Engine (GEE) for 2014, 2016, and 2017.For the remote sensing analysis, surface reflectance Landsat 8 imagery was used for the study area (path 19, row 49).Orthorectification and coregistration were completed on all the imagery used in this study by USGS prior to integration in the GEE (Google, Inc., MountainView, CA) platform.The imagery was atmospherically corrected using LaSRC software.Clouds and cloud shadows were masked with the CFmask function.Image composites were created by selecting imagery with the lowest possible cloud cover (<10%) using the composite algorithm available in GEE.False color composites and normalized difference vegetation indexes (NDVI = (NIR − RED)/(NIR + RED)) were computed to enhance the interpretability of the imagery.
A stratified random sampling design was used to select training and validation data to perform the supervised classification of Landsat imagery [34].Of the 1000 sampling points, 700 were randomly selected as training data and the remaining 300 were used as validation points.The training and validation data were primarily created from field surveys that were conducted within the MGL to collect ground-truthed GPS points of land cover types.Training points located in remote areas that could not be visited in the field due to safety or cost constraints were visually interpreted using freely The MGL is becoming increasingly threatened as unsustainable land use practices reduce the land's capacity to provide life-supporting ecosystem services.The region is farmed predominantly through slash-and-burn agriculture.Traditionally, farmers will cultivate a plot until it decreases in productivity, at which point it will be left to re-grow natural vegetation for a period of ten to fifteen years.During this fallow period, the soil is able to regain fertility for subsequent cultivation.However, in the last few decades, the fallow period of most plots has been reduced to two or three years due to a shortage of land brought about by a steady increase in population [7,8].Consequently, the soil is not able to completely regain its fertility, resulting in more numerous and shorter agricultural cycles and increased deforestation [9].While farmers continue to clear secondary-growth forests that were left fallow, they have also begun to cultivate mature forests that have not been cleared in the recent past.There has also been an increase in other unsustainable agricultural practices in the MGL such as large-scale citrus and banana plantations [10].
Protected area and sustainable livelihood managers in Belize and throughout the tropics have limited resources.Therefore, detailed information on the location of sites vulnerable to forest conversion is necessary to prioritize areas for law enforcement and compliance, sustainable management, and community outreach.However, the location and size of sites most vulnerable to future tropical deforestation remain unknown.Thus, a need exists to identify the potential future spatial distribution and magnitude of forest loss to strategically implement conservation efforts where they will be the most effective.
Data-driven models of deforestation patterns can assist forest protection and management organizations in conservation planning to efficiently allocate resources and produce the greatest conservation impact.A plethora of research has focused on the locations and rates of forest cover change based on remote sensing technology.Recently, the results of these studies have been used to identify the predictors of change and to assess specific areas where forest loss is likely to occur in the future [11][12][13][14][15][16][17][18][19][20].
Deforestation vulnerability and future forest cover change can be predicted using simulation models and empirical models [21].These models assess how and why changes occurred in the past to understand the main drivers of deforestation in the present and therefore predict where and how much change will arise in the future.The most common methods used in land cover prediction modeling are weights of evidence [18], logistic regression [14,22], and multi-layer perceptron (MLP) neural networks [23].Recent research has shown that MLP neural networks perform better than other methods in land change modeling since they can more adequately process non-linear relationships [11,24,25].The algorithm for the MLP neural network is non-parametric and does not consider multicollinearity.In addition, it allows for the fitting of complex data structures even if insignificant variables are included in the model [17,26].This study employed a MLP neural network model to determine the sites most susceptible to future deforestation based on the anthropogenic conversion of forest in the MGL using remote sensing and spatial variables that represent the regional drivers of deforestation.In this research, vulnerable areas are defined as regions that are susceptible to forest loss as a result of human activities, primarily agriculture.
Landscape-scale forest cover change studies have been previously published for Belize [27][28][29], including the southern state of Toledo [30,31], but the data are six years old or more.Additionally, although models have been developed to predict the potential spatial drivers of deforestation in Belize [29,32], no assessment has been published that predicts the magnitude, spatial extent and patterns of forest loss for any portion of the country.Fine scale assessments are needed by protected area and livelihood managers at the regional landscape level since the drivers of deforestation are typically fine scale and site-specific [33].This method could be applied more widely to conservation planning in other landscapes, especially regions with severe financial constraints, as a way to prioritize the protection of threatened areas.

Forest Cover Classification and Change Analysis
Forest cover was classified utilizing a supervised classification of Landsat 8 Operational Land Imager (OLI) data in Google Earth Engine (GEE) for 2014, 2016, and 2017.For the remote sensing analysis, surface reflectance Landsat 8 imagery was used for the study area (path 19, row 49).Orthorectification and coregistration were completed on all the imagery used in this study by USGS prior to integration in the GEE (Google, Inc., MountainView, CA) platform.The imagery was atmospherically corrected using LaSRC software.Clouds and cloud shadows were masked with the CFmask function.Image composites were created by selecting imagery with the lowest possible cloud cover (<10%) using the composite algorithm available in GEE.False color composites and normalized difference vegetation indexes (NDVI = (NIR − RED)/(NIR + RED)) were computed to enhance the interpretability of the imagery.
A stratified random sampling design was used to select training and validation data to perform the supervised classification of Landsat imagery [34].Of the 1000 sampling points, 700 were randomly selected as training data and the remaining 300 were used as validation points.The training and validation data were primarily created from field surveys that were conducted within the MGL to collect ground-truthed GPS points of land cover types.Training points located in remote areas that could not be visited in the field due to safety or cost constraints were visually interpreted using freely available high-resolution satellite imagery and aerial imagery from Google Earth TM .Ecosystem layers [35] and fire point data [36] were used as reference materials to help interpretation, especially of non-forest natural areas such as savannas (including transitional zones), wetlands (including mangroves), and large bodies of water.
Three classes of land cover were classified: forests, anthropogenic areas, and non-forest natural areas.Anthropogenic areas include small-to-medium scale slash-and-burn agriculture, villages or settlements, large-scale citrus and banana plantations, cattle pastures, and medium-scale Mennonite community agriculture.To model and predict the three land cover classes, a Classification and Regression Trees (CART) classifier was deployed in GEE [37].The land cover classification data from 2014 to 2016 were used to perform a forest cover change analysis.Regions classified as non-forest natural areas were excluded from the forest cover change analysis, since the focus of the study was on the conversion of broadleaf tropical forest to anthropogenic areas.Additionally, Landsat 8 imagery was classified for 2017 in GEE to conduct the validation of the prediction model.
The extremely dynamic and highly successional nature of the patches that are continually subjected to slash-and-burn agriculture in the MGL present a unique challenge to determining forest cover change.Within this patchy landscape, regenerating secondary-growth forest only lies in fallow for two to three years before it is subsequently converted to agriculture once again.Thus, calculating such highly dynamic rates of deforestation and natural regeneration is generally irrelevant based on the overall balance throughout this landscape.In this analysis, early and mid-regeneration forests have been excluded from the calculation of deforestation within the MGL.Therefore, the deforestation calculated for this study only represents deforestation that has occurred in mature forest.This follows the methodology of other forest cover change analysis in Belize in which the researchers defined "forest" as natural mature forests and did not include secondary-growth forests in their studies [27,28].
To calculate deforestation for only mature forests, the forest class that was described from the Landsat imagery was subdivided into early-and-mid-regeneration forests (i.e., areas that have been continually subjected to slash-and-burn agriculture since 1980) and mature forests.Early-andmid-regeneration forests were distinguished from mature forests based on the supervised classification of Landsat MSS, TM, ETM+, and OLI imagery from 1980 to 2013 in ENVI and GEE from previous studies of unpublished data [38].In these prior studies, anthropogenic areas were classified based on Landsat imagery from 1980 to 2013 (Figure 2).All regions that were classified as historical anthropogenic areas in the previous study but were classified as forest in this study based on the 2014 Landsat imagery were re-classified as early-and-mid-regeneration forests for the 2014 land cover classification.The remaining forest cover class areas were re-classified as mature forest.The same classification scheme was replicated to create early-and-mid-regeneration forests for 2016 and 2017.
To evaluate accuracy of land cover classification, a confusion matrix [39] was generated using 300 validation points.These data were collected from ground-truthed GPS points as well as freely available high-resolution satellite imagery and aerial imagery from Google Earth TM .Producer's accuracy (the probability that a value in a class was classified correctly) was calculated by the total number of correctly classified pixels divided by the total number of pixels truly in that class.User's accuracy (the probability that a value predicted to be in a certain class truly is that class) was calculated by the number of correctly classified pixels divided by the total number of pixels predicted within that class.The overall accuracy was also calculated by summing the number of correctly classified values and dividing by the total number of values.

Forest Cover Change Prediction Model
To assess vulnerability to future deforestation and predict the forest patches that may be cleared within the next ten years, a tool called the Land Change Modeler (LCM) [40] was implemented from within the TerrSet Geospatial Monitoring and Modeling Software.The LCM analyzes previous land cover change and spatial variables that represent regional drivers of deforestation, models the deforestation vulnerability, and predicts the patches that may be cleared in the future.The model computes the following steps: (i) historic change analysis; (ii) spatial deforestation driver suitability evaluation and selection; (iii) transition potential map creation; (iv) land demand and allocation estimation; (v) validation; and (vi) future land cover vulnerability and prediction.The process is described below.

Historic Change Analysis
In Step 1, another forest cover change analysis from 2014 to 2016 was conducted in LCM to identify focal areas of change in LCM.This was conducted in addition to the forest cover change analysis in GEE in order to record the historic patterns of change in LCM.

Spatial Deforestation Driver Suitability Evaluation and Selection
In Step 2, spatial drivers of deforestation were evaluated and selected to input into the model for the creation of the transition potential maps.Within the MGL, the clear-cutting of trees for timber in broadleaf forest is currently rare and the main driver of deforestation is agriculture Farmers in Belize are more likely to clear a patch of mature forest if it is easily accessible (near a road or settlement), easily cleared (located near the edge of the forest), suitable for agriculture, and not located in a strict protected area [29,32].Previous studies in which models identified the major drivers of deforestation in tropical forests for Latin America found similar results to models conducted for Belize.For example, higher deforestation occurred in close proximity to roads, settlements, forest

Forest Cover Change Prediction Model
To assess vulnerability to future deforestation and predict the forest patches that may be cleared within the next ten years, a tool called the Land Change Modeler (LCM) [40] was implemented from within the TerrSet Geospatial Monitoring and Modeling Software.The LCM analyzes previous land cover change and spatial variables that represent regional drivers of deforestation, models the deforestation vulnerability, and predicts the patches that may be cleared in the future.The model computes the following steps: (i) historic change analysis; (ii) spatial deforestation driver suitability evaluation and selection; (iii) transition potential map creation; (iv) land demand and allocation estimation; (v) validation; and (vi) future land cover vulnerability and prediction.The process is described below.

Historic Change Analysis
In Step 1, another forest cover change analysis from 2014 to 2016 was conducted in LCM to identify focal areas of change in LCM.This was conducted in addition to the forest cover change analysis in GEE in order to record the historic patterns of change in LCM.

Spatial Deforestation Driver Suitability Evaluation and Selection
In Step 2, spatial drivers of deforestation were evaluated and selected to input into the model for the creation of the transition potential maps.Within the MGL, the clear-cutting of trees for timber in broadleaf forest is currently rare and the main driver of deforestation is agriculture Farmers in Belize are more likely to clear a patch of mature forest if it is easily accessible (near a road or settlement), easily cleared (located near the edge of the forest), suitable for agriculture, and not located in a strict protected area [29,32].Previous studies in which models identified the major drivers of deforestation in tropical forests for Latin America found similar results to models conducted for Belize.For example, higher deforestation occurred in close proximity to roads, settlements, forest edges, and areas with a low level of protection in many previous studies [41][42][43][44][45].In addition, forests with fertile soils were more likely to be converted to agriculture than forests containing nutrient poor soils [46,47].Therefore, the following deforestation drivers were selected as spatial variables in this study: proximity to roads, proximity to settlements, proximity to forest edges, agricultural suitability, and level of protection.Other factors that influenced the selection of spatial variables were spatial data availability and the visual inspection of forest cover change in the MGL since 1980.The MLP neural network mitigates the possibility of correlated variables and data redundancy and therefore these issues are minor in this study [48].
The relationship between the spatial variables used in this study and predicted forest changes were tested using Cramer's V coefficient to evaluate their suitability for inclusion in the model.Eastman [40] suggests that all deforestation drivers with a Cramer's V coefficient less than 0.15 should be discarded.All spatial drivers that were included in the model in this study had a strong Cramer's V predictive power (V ≥ 0.3) and significant p value (p < 0.05).Biophysical characteristics such as slope, elevation, and soil type were tested in the LCM as potential deforestation drivers but were determined to not be strong predictors to change.Therefore, these biophysical datasets were not selected as spatial drivers of deforestation for this model.
The raster layers for the spatial variables of settlement edge and forest edge were created from the land cover classification rasters in this study.The raw road dataset was developed by Meerman and Clabaugh [35].All proximity-based raster data were created from the aforementioned raw data, utilizing Euclidean distance models in the LCM.
The agricultural suitability dataset was created by King et al. [49], based upon field surveys of soil fertility, soil type, flood risk, erosion potential, salinity, geology, and slope.Although the agricultural suitability index was calculated based on field surveys from over thirty years ago, the spatial data were examined by agricultural extension officers and found to be accurate and adequate based on their extensive local knowledge.The raw protected area dataset was developed by Meerman and Clabaugh [35].The variables that represent the level of protection were ranked loosely based on the International Union for Conservation of Nature's (IUCN's) protected area categories, which are classified according to management goals and the level of permitted human activity.For example, the following are designated levels of protection from least likely to most likely to be deforested: (i) strict nature reserves; (ii) national parks and wildlife sanctuaries; (iii) agroforestry concessions; (iv) forest reserves; and (v) areas with no protection.All of the spatial data variables were manipulated in ArcGIS 10.5 [50] and LCM prior to use in the model.The spatial variable rasters were created with a cell size of 30 m to match the resolution of the land classification rasters.

Transition Potential Map Creation
In Step 3, transition potential maps were created to represent the likelihood for a patch of mature forest to be converted to an anthropogenic patch [24].These were generated using data from the forest cover change analysis and spatial variables that had been selected in step two as the drivers of deforestation in the MGL.Transition potentials can be developed with several methods including logistic regression [14,22], multi-layer perceptron (MLP) neural networks [23], and weights of evidence [18].This model was implemented with a MLP neural network, as recent studies have found that it outperforms other methods [11,24,25].This study achieved an accuracy training rate of 82%, which is above the minimum required level.A detailed description of the training procedure is documented by [24,40].

Land Demand and Allocation Estimation
The LCM predicts forest change allocation based on a multi-objective land allocation (MOLA) algorithm [51] in Step 4 of the process.This study applied a Markov chain model in LCM to estimate the quantity of change that could develop at a specified year in the future based on previous land cover data and the probability of conversion from one land cover to another [40].A prediction map for 2017 was generated in LCM using the Markov chain model to validate the model's accuracy.

Model Validation
The predicted output of a model should be compared to observed data to validate and assess the model's predictive power.To validate this model in Step 5, the area under the curve (AUC) of the Receiver Operating Characteristic (ROC) was calculated to evaluate the goodness of fit, following the authors of [14,17,22,24,46,52].This technique involves a series of comparisons calculated in the Validate Module within the LCM software.The procedure compares the reference forest cover change map with the model output map and evaluates the proportion of "false alarms" and "hits" [24].This generates the AUC figure, with a value of one representing a perfect model and a value of 0.5 representing an accuracy equivalent to a random fit.The 2017 prediction map generated in Step 4 was used to assess the accuracy of the forest cover change prediction model, since it could be compared to the 2017 land cover classification map (Figure 3) and the forest cover change analysis that was conducted in this study.

Deforestation Vulnerability and Forest Cover Change Prediction
After the model was validated, Step 6 was performed, producing a map of the vulnerability to future forest conversion in the MGL with values that represent low (0) to high ( 1) vulnerability.The model also produced a prediction map for the land cover scenario of 2026.The vulnerability map identifies all of the locations for the entire landscape that are susceptible to deforestation at any point in time and their level of vulnerability, whereas the prediction map determines the most likely land cover classifications for a specific year.To predict the magnitude, spatial extent, and patterns of forest loss for 2026, the original land cover map of 2016 was intersected with prediction map for 2026 and the projected deforestation for the next ten years was calculated in ArcGIS.

Land Cover Classification and Model Validation
The land cover classification achieved an overall accuracy of 88% for 2014, 94% for 2016, and 95% for 2017.Class specific producer's and user's accuracies ranged from 86% to 96%.This exceeds thresholds previously established as acceptable, e.g., the USGS suggests a threshold of 85% [53] and a threshold of 80% is suggested in [54].This prediction model achieved a goodness of fit of 84% from the AUC of the ROC.Thus, the model is a viable predictor of change.

Forest Cover Change Analysis
The land cover classification analysis for 2016 shows that the MGL is currently a highly forested, intact landscape (Figure 4).As of 2016, a total of 75% of the MGL remained in mature forest.Over 90% of the MGL currently retains natural vegetation (including savanna, wetland, and early to mid-regeneration forests in fallow), which was calculated with the data from this analysis and data from [35].However, the forest cover change analysis shows that this landscape is threatened by deforestation.From 2014 to 2016, 2090 ha (5165 ac) of mature forest were cleared.

Forest Cover Change Prediction Model
The maps depicting the vulnerability to future deforestation, the predicted forest cover for 2026, and predicted mature forest conversion for 2026 are presented in Figures 5-7.The assessment indicates that the agricultural frontier will continue to expand into mature forests.According to the model, about 7392 ha (18,267 ac) of mature forest are expected to transition to anthropogenic areas, such as slash-and-burn agriculture and plantations (Figure 6).It predicts that mature forests will decrease from 75.0% forest cover in 2016 to 71.9% in 2026 for the entire MGL.No deforestation is projected for any of the reserves with a strict level of protection (GSCP, BNR, CBWS, and PCNP); while a slight increase in forest cover loss is projected for the forest reserves (CRFR, MMNFR, DRFR, and SBFR).The forest within the agroforestry concession located in MMNFR is projected to naturally regenerate.The majority of the deforestation is predicted to occur in the community zone, with only slightly over 28,730 ha (71,000 ac) of mature forest predicted to remain in this region by 2026 (Figures 6-7).While the mature forest cover in protected areas will likely remain around 84.0% over the next ten years, mature forests in the community zone will likely decrease from 49.5% to 40.0%, according to the model.This is also predicted to result in changes of spatial distribution of mature forests in the community zone, such as increased fragmentation.The main regions of future predicted deforestation are the de-reserved zones of MMNFR, which are likely to greatly decrease in the next few years.The model predicts that the forest cover will decrease to 29.2% in the 2006 de-reserved zone and to 27.8% in the 2015 de-reserved zone.In terms of structural connectivity, the MGL is dominated by core intact mature forest, most of which is primarily located in protected areas.Several protected areas in the MGL, such as Bladen Nature Reserve (BNR), Golden Stream Corridor Preserve (GSCP), Cockscomb Basin Wildlife Sanctuary (CBWS) and Payne's Creek National Park (PCNP), did not experience deforestation from 2014 to 2016, while others did, such as Columbia River Forest Reserve (CRFR), Maya Mountain North Forest Reserve (MMNFR), Deep River Forest Reserve (DRFR), and Swasey-Bladen Forest Reserve (SBFR), but usually at the edges of their borders.Most of the mature forest cover loss that has occurred within the MGL is located outside of the region's reserves, referred to in this study as the community zone.From 2016 to 2017, the annual deforestation rate within protected areas in the MGL was only 0.04%, while the deforestation rate in the community zone of the MGL was 0.90%.
Two of the main regions of deforestation are the zones that have been de-reserved from MMNFR.About 1315 ha (3250 ac) of MMNFR were de-reserved in 2006 (referred to as the 2006 de-reserved zone) adjacent to its current southern boundary, and another 935 ha (2310 ac) were de-reserved in 2015 (referred to as the 2015 de-reserved zone), adjacent to its current southeastern corner.The current mature forest cover of the 2015 de-reserved zone is 63.9% and it was approximately 76.2% prior to the elimination of its protection.The annual deforestation rates in 2016 (4.35%) and 2017 (4.54%) were more than double the rate from the year immediately prior to downsizing (1.71%).The current mature forest cover of the 2006 de-reserved zone is only 50%.While additional analyses would need to be conducted to determine the forest cover and deforestation rates prior to downsizing, this figure is much lower than the current mature forest cover of 97.7% within the entire MMNFR.

Forest Cover Change Prediction Model
The maps depicting the vulnerability to future deforestation, the predicted forest cover for 2026, and predicted mature forest conversion for 2026 are presented in Figures 5-7.The assessment indicates that the agricultural frontier will continue to expand into mature forests.According to the model, about 7392 ha (18,267 ac) of mature forest are expected to transition to anthropogenic areas, such as slash-and-burn agriculture and plantations (Figure 6).It predicts that mature forests will decrease from 75.0% forest cover in 2016 to 71.9% in 2026 for the entire MGL.No deforestation is projected for any of the reserves with a strict level of protection (GSCP, BNR, CBWS, and PCNP); while a slight increase in forest cover loss is projected for the forest reserves (CRFR, MMNFR, DRFR, and SBFR).The forest within the agroforestry concession located in MMNFR is projected to naturally regenerate.The majority of the deforestation is predicted to occur in the community zone, with only slightly over 28,730 ha (71,000 ac) of mature forest predicted to remain in this region by 2026 (Figures 6 and 7).While the mature forest cover in protected areas will likely remain around 84.0% over the next ten years, mature forests in the community zone will likely decrease from 49.5% to 40.0%, according to the model.This is also predicted to result in changes of spatial distribution of mature forests in the community zone, such as increased fragmentation.The main regions of future predicted deforestation are the de-reserved zones of MMNFR, which are likely to greatly decrease in the next few years.The model predicts that the forest cover will decrease to 29.2% in the 2006 de-reserved zone and to 27.8% in the 2015 de-reserved zone.

Community Zone
Although the MGL is still mostly forested, the agricultural frontier continues to advance into mature forests due to unsustainable small-scale and large-scale agriculture.Farmers have implemented shorter crop cycles, which has led to a decrease in soil fertility and the expansion of agriculture land into mature forest.The model predicts that this will continue to occur in the next ten years.The decreased area of remaining mature forest projected in the community zone will most likely increase pressure on the farmers' ability to find fertile agricultural plots in this region and may result in increased pressure on reserves beyond 2026.
The zones that have been de-reserved from MMNFR are most likely losing forest cover rapidly due to a lack of proper landscape management.The de-reservation of the region is an event called downsizing, in which an official decree legally changes a protected area's boundaries, resulting in a reduction in size.Although reserves are assumed to be permanent features within the landscape, research suggests that protected area downgrading, downsizing, and degazettement (PADDD) is actually common and widespread and can result in adverse conservation impacts such as forest cover loss [55,56].This model was created under the assumption that protected area boundaries would not change in the next ten years.If any reserves are downgraded, downsized, or degazetted, these areas may follow the same trajectory as the MMNFR de-reserved zones and would likely be more vulnerable to future deforestation than the model predicted.

Protected Areas
The status and remoteness of protected areas most likely plays an important role in the level of impact projected within their boundaries.Although no deforestation is predicted to occur at GSCP, a private protected area, the edges are highly vulnerable to future forest conversion.This is most likely because it is located directly adjacent to communities and agricultural plots and a road bisects the property.A very slight increase in deforestation is projected for Boden Creek Ecological Preserve (BCEP), a private protected area, most probably because it is adjacent to drivers of deforestation such as roads, settlements, and forest edges.No deforestation is predicted to occur in BNR and its boundaries are not vulnerable to deforestation, most likely due to its remote location and its strict regulations in which research, education and law enforcement are the only permitted human activities.The projected increase in deforestation along the boundaries of all four forest reserves, CRFR, MMNFR, DRFR, and SBFR, may be related to their protection status and close proximity to the drivers of deforestation.
Forest cover loss has been concentrated in CRFR for many years in two different regions: the southern and western edges.On the western boundary of CRFR, which lies adjacent to the Guatemala border, small clearings have advanced into the reserve, and the model predicts that this will most likely progress in the future.The Guatemalan side of the international border near CRFR has been heavily deforested [57].Guatemalan citizens began crossing the border into Belize in the early 1990s to utilize the land for illegal farming, cattle ranching, marijuana production, timber and non-timber forest product harvest, and game hunting [58][59][60].The remoteness of the border, lack of personnel, lack of financial resources, and high danger of armed encounters are barriers to enforcement by the Belize Defense Force and Forest Department [60].In a different section located near the southern boundary of CRFR, an old logging road has provided access to an agricultural area within the reserve that has expanded over time, which will continue in the future according to the model.In addition to the expansion of current cleared areas, the model also predicts new incursions, especially along the southern boundary of CRFR.Without proper conservation planning and strategic placement of patrols on the southern boundary, these forests could be lost.
Deforestation has occurred in the southeast corner of MMNFR for several years.This section is currently a cacao agroforestry concession but, prior to 2015, the boundary contained a logging concession.Clear-cutting was practiced in this site prior to the change in management of the concession, which is reflected in the deforestation exhibited in this area prior to 2015.The forest has already begun to regenerate and will be allowed to continue to regrow after the cacao reaches full production in approximately four years.The model was implemented with the consideration that the level of permitted human activity in an agroforestry concession would most likely have a lower impact on intact forest cover than a logging concession.It predicts that no mature forest cover loss will occur in the next ten years under this scenario and the clearings will actually naturally regenerate.Although the model does predict that this area is vulnerable to future deforestation, this is most likely influenced by its past land use history and its close proximity to drivers of deforestation such as roads and agricultural plot edges in the 2015 de-reserved zone.

Limitations
The model of deforestation vulnerability and predicted deforestation extent, magnitude, and spatial distribution is based on several assumptions.Vulnerability and future forest patterns were based on historical deforestation patterns and it is assumed that the forest will change in the same manner as the forest cover change analysis of 2014-2016.While recent historical patterns can be fairly accurate predictors of change, tropical deforestation is a dynamic process and uncertainties can affect the trajectory of tropical forests.It does not account for climate change or possible future natural disasters such as hurricanes, which could cause massive deforestation.The model was also implemented under the assumption that the drivers of deforestation will not vary significantly within the next ten years.It implies that no additional roads will be constructed and no protected areas will be established or de-reserved.This study should be recognized as a baseline model that can be monitored and adapted by iterative assessments.The model was designed to be updated with new data should a new situation or driver arise.For instance, the model could be altered to assess the impact of the de-reservation of a protected area or the building of a new road.This assessment was also limited by the availability of spatial data.Some data that were significant drivers of deforestation in other studies [29,43], such as population density and climatic factors (annual rainfall and dry-season intensity), could not be tested as potential drivers of deforestation since these data are not available in a spatial form for southern Belize.This model could be modified if new or updated spatial datasets become available in the future.

Implications for Informing Conservation Planning
The results of this study have helped to inform the conservation planning process of stakeholders within the MGL.The vulnerability and prediction maps are currently being used by protected area and sustainable livelihood managers to identify and prioritize where to strategically focus conservation actions in an effort to prevent future deforestation and alternatively retain an intact, highly forested landscape.The locations found to be the most vulnerable to forest conversion are priority locations for conservation actions such as encouraging the adoption of sustainable agroforestry practices, implementing community outreach and increasing protection efforts.Law enforcement and compliance actions, such as increased patrols, will be implemented within the most threatened regions of protected areas.Additionally, community outreach and sustainable agricultural practices will be implemented in the communities that are the most vulnerable to deforestation to prevent future forest conversion.Trainings on proper fire management will be given to slash-and-burn agricultural farmers in regions that are the most vulnerable to forest conversion in an effort to reduce the risk of escaped fires.Outreach on agroforestry and inga alley cropping will be given to farmers in communities that are the most threatened by forest loss.By shifting from slash-and-burn agriculture to agroforestry or inga alley cropping, farmers can increase the soil fertility on their land and reduce the tendency to cut mature forest.

Conclusions
This study employed a multi-layer perceptron (MLP) neural network model to predict the sites most susceptible to future deforestation based on recent forest cover change in the Maya Golden Landscape (MGL) and spatial variables that represent the regional drivers of deforestation.The MGL is unique in that it has remained a highly forested, intact landscape with 75% of its cover in mature forest.However, challenges exist in retaining this standing forest.
The model predicts that the MGL will continue to exhibit an expansion of agriculture into mature forest, with a decrease of mature forest cover to 71.9% in 2026.This is most likely due to the implementation of shorter crop cycles which have decreased soil fertility.Most of the mature forest cover is projected to be lost in the community zone.The areas that have been de-reserved from MMNFR are some of the most vulnerable regions with a projected future forest cover of less than 30%, most likely a result of unsustainable landscape management.The model predicts that no deforestation will occur within most strict protected areas, most likely due to their remote location and protection status.A slight increase in forest cover loss is projected to occur along the boundaries of all four forest reserves, most likely associated with their close proximity to the drivers of deforestation and their level of permitted human activity.Although only a slight increase in deforestation is projected to occur in reserves, the potential high level of deforestation in the community zone may put pressure on the protected areas beyond 2026.The vulnerability and prediction maps are currently being used to strategically focus conservation efforts by stakeholders to effectively allocate resources in an attempt to retain an intact, highly forested landscape.

Figure 1 .
Figure 1.Location of study area and protected areas of the Maya Golden Landscape.

Figure 1 .
Figure 1.Location of study area and protected areas of the Maya Golden Landscape.

Figure 3 .
Figure 3. Land cover of the Maya Golden Landscape in 2017: (a) as classified using Landsat 8 in Google Earth Engine; and (b) as predicted by the Land Change Modeler.Non-forest natural areas include savannas, wetlands, large bodies of water, and mangroves.

Figure 3 .
Figure 3. Land cover of the Maya Golden Landscape in 2017: (a) as classified using Landsat 8 in Google Earth Engine; and (b) as predicted by the Land Change Modeler.Non-forest natural areas include savannas, wetlands, large bodies of water, and mangroves.

Figure 4 .
Figure 4. Land cover classification of the Maya Golden Landscape in 2016.Non-forest natural areas include savannas, wetlands, large bodies of water, and mangroves.

Figure 4 .
Figure 4. Land cover classification of the Maya Golden Landscape in 2016.Non-forest natural areas include savannas, wetlands, large bodies of water, and mangroves.

Figure 5 .
Figure 5. Vulnerability to future deforestation, higher values represent higher vulnerability.Figure 5. Vulnerability to future deforestation, higher values represent higher vulnerability.

Figure 5 .
Figure 5. Vulnerability to future deforestation, higher values represent higher vulnerability.Figure 5. Vulnerability to future deforestation, higher values represent higher vulnerability.

Figure 5 .
Figure 5. Vulnerability to future deforestation, higher values represent higher vulnerability.

Figure 6 .
Figure 6.Forest cover and anthropogenic areas predicted for 2026.Figure 6. Forest cover and anthropogenic areas predicted for 2026.