Modelling Global Deforestation Using Spherical Geographic Automata Approach

: Deforestation as a land-cover change process is linked to several environmental problems including desertiﬁcation, biodiversity loss, and ultimately climate change. Understanding the land-cover change process and its relation to human–environment interactions is important for supporting spatial decisions and policy making at the global level. However, current geosimulation model applications mainly focus on characterizing urbanization and agriculture expansion. Existing modelling approaches are also unsuitable for simulating land-cover change processes covering large spatial extents. Thus, the objective of this research is to develop and implement a spherical geographic automata model to simulate deforestation at the global level under different scenarios designed to represent diverse future conditions. Simulation results from the deforestation model indicate the global forest size would decrease by 10.5% under the “business-as-usual” scenario through 2100. The global forest extent would also decline by 15.3% under the accelerated deforestation scenario and 3.7% under the sustainable deforestation scenario by the end of the 21st century. The obtained simulation outputs also revealed the rate of deforestation in protected areas to be considerably lower than the overall forest-cover change rate under all scenarios. The proposed model can be utilized by stakeholders to examine forest conservation programs and support sustainable policy making and implementation.


Introduction
Deforestation is a land-cover change (LCC) process caused by natural and anthropogenic factors further entailing environmental degradation with several negative consequences both at the regional and global scales [1,2]. Rates of deforestation especially in developing countries, triggered by factors such as agricultural expansion, timber production, forest fire, mining, and urbanization have been increasing over the last century [3][4][5]. These trends of decreasing forest cover and deteriorating conditions have resulted in deforestation becoming a major global environmental issue considering the several critically important ecosystem services and functions forests provide [6,7]. Forests are important areas for biodiversity, with approximately 80% of the world's terrestrial biodiversity found in forest regions [8]. Further, forests represent the largest terrestrial sink of carbon dioxide (CO 2 ) and are globally responsible for significant carbon stocks [9,10]. These highlight the important roles forests play in the Earth's biochemical and ecological systems.
Geosimulation modelling has become an important tool for representing land-cover change (LCC) processes and assisting in understanding the interaction between anthropogenic activities and the impact of deforestation on other environmental systems, permitting spatial analyses of the underlying causes of this dynamic process [11]. Further, the simulation of possible LCC scenarios provides a useful mechanism to inform environmental and forest management policies and decision making for providing valuable insights for ISPRS Int. J. Geo-Inf. 2023, 12, 306 2 of 21 developing appropriate measures to alleviate the negative impacts of deforestation [12]. Specifically, data on forest cover and future trajectories provide significant information for estimating carbon stocks, ecosystem service evaluation, and forestry conservation [13,14]. Assessing the performance and effectiveness of environmental policies such as the Reduction in Emissions from Deforestation and Forest Degradation (REDD+) policies requires detailed spatial data on forest-cover change and future scenarios [15].
Forests can be considered as complex biophysical spatial systems with many components, and some considerably depend on human interactions at the local level to give rise to global patterns of deforestation processes over time [16]. Thus, geosimulation modelling approaches are seen as suitable for representing forest change processes. Accordingly, several geosimulation models have been implemented to represent the dynamics of forest changes as a complex spatial process including approaches based on cellular automata (CA) [17,18], and some were enhanced with techniques such as Markov chain [19,20], logistic regression [21,22], multi-criteria evaluation (MCE) [23,24], machine learning [25], and deep learning [26]. Several studies have also been incorporating human interactions to represent deforestation processes using agent-based geosimulation models (ABMs) [27,28]. However, these modelling approaches are all developed mainly to operate on small spatial extents and implemented to simulate forest-cover change dynamics at the local and regional levels.
The use of existing geosimulation models at larger spatial scales presents challenges that are peculiar to spatial modelling at the global level. Primarily, these models do not consider the curvature of the Earth's surface when modelling at larger extents, which can lead to errors in spatial and statistical analyses due to spatial distortions caused by planar map projections [29]. The limitations of using planar spatial models for analyses and simulations at the global level have been documented in the scientific literature [30][31][32], with spherical models proposed as a possible solution. While spatially explicit land-cover change models have become prevalent over the last decade, geosimulation models for deforestation are still scarce, with existing global applications typically focusing on simulating urbanization [33,34] and agricultural expansion [35]. In order to improve simulation results, CA models are often integrated with other techniques such as spatial multi-criteria evaluation (MCE) to identify suitable or susceptible locations for the potential occurrence of the geographic phenomena and then guide transition rules. The MCE technique provides a comprehensive approach that combines several often-conflicting criteria based on suitability functions, weights, and their overall aggregation [36]. The approach has been implemented in several applications, including land-cover change [37,38], deforestation [23], and urban growth [39], although all these studies address small spatial extents. Therefore, the main objective of this research study is to develop and implement a spherical geographic automata (SGA) modelling approach by integrating MCE and cellular automata to simulate the process of deforestation at the global level and considering the curved surface of the Earth.

Spherical Deforestation Model Overview
The methodology extends the theoretical concepts of the spherical geographic automata (SGA) approach [40] and integrates susceptibility analysis for global deforestation modelling. The methodological flow chart of the modelling approach is presented in Figure 1. The spherical component of the model is operationalized with the use of a discrete global grid system (DGGS) [41] and hexagonal spatial tessellations as a base unit that allows for geospatial data representation at the global level and with consideration of the curvature of the Earth's surface. The MCE technique is used to identify susceptible locations for forest-cover loss using several criteria as possible drivers of deforestation. Three scenarios have been developed to represent possible future deforestation processes under different conditions.

Global Deforestation Spherical Geographic Automata
The spherical geographic automata (SGA) component is the central part of the proposed modelling methodology, and it is designed to simulate the process of deforestation at the global level. The SGA utilizes a spherical cell space based on DGGS and comprising hexagonal tessellation covering the Earth's curved surface. As a geospatial model, DGGS applies a spherical grid framework to partition and represent the curvature of the Earth's surface [41]. The DGGS spatial model is based on an icosahedron polyhedron with equalarea hexagonal cells. When used to tesselate spherical surfaces, hexagonal cells are the most compact and offer uniform adjacency and neighbouring relationships over other regular polygons such as squares and triangles. The global spatial datasets are transformed into hexagonal spatial tessellations as the model input. The spherical geographic deforestation model extends the previous research study [42] and can be formulated as follows: where ℎ +1 is the state of the hexagonal cell h at the next time step t+1, ℎ denotes the state of the hexagonal cell at initial time t, ℎ represents the hexagonal neighbourhood of six cells surrounding the central cell, def_ is the overall susceptibility value obtained for each hexagonal cell, is the function of transition rules that determine how the state of cells changes over time, and Δ is the discrete time step representing one iteration of the model. The effects of protected areas on deforestation as constraints and values of susceptibility analysis of deforestation are also considered. The function of transition rules represents the actual dynamics of the deforestation process. During each iteration,

Global Deforestation Spherical Geographic Automata
The spherical geographic automata (SGA) component is the central part of the proposed modelling methodology, and it is designed to simulate the process of deforestation at the global level. The SGA utilizes a spherical cell space based on DGGS and comprising hexagonal tessellation covering the Earth's curved surface. As a geospatial model, DGGS applies a spherical grid framework to partition and represent the curvature of the Earth's surface [41]. The DGGS spatial model is based on an icosahedron polyhedron with equalarea hexagonal cells. When used to tesselate spherical surfaces, hexagonal cells are the most compact and offer uniform adjacency and neighbouring relationships over other regular polygons such as squares and triangles. The global spatial datasets are transformed into hexagonal spatial tessellations as the model input. The spherical geographic deforestation model extends the previous research study [42] and can be formulated as follows: where GA t+1 h is the state of the hexagonal cell h at the next time step t + 1, GA t h denotes the state of the hexagonal cell at initial time t, HN t h represents the hexagonal neighbourhood of six cells surrounding the central cell, S t def_global is the overall susceptibility value obtained for each hexagonal cell, f is the function of transition rules that determine how the state of cells changes over time, and ∆T is the discrete time step representing one iteration of the model. The effects of protected areas on deforestation as constraints and values of susceptibility analysis of deforestation are also considered. The function of transition rules represents the actual dynamics of the deforestation process. During each iteration, cells representing forest are converted to the dominant non-forest land-cover type based on the cell's neighbourhood, susceptibility value, and constraint parameter.

Datasets
The study area in this research study encompasses the entire global land surface except for Antarctica, and several spatial datasets with global extent were acquired to implement the model. Land-cover datasets were obtained from the European Space Agency (ESA) portal [43], global roads dataset from the Global Roads Inventory Project (GRIP) [44], protected areas from the World Database on Protected Areas (WDPA) [45], past forest disturbance from the Global Wildfire Information System (GWIS) [46], population density from the LandScan portal [47], and elevation dataset from the United States Geological Survey (USGS) portal [48]. All spatial datasets were converted into Icosahedral Snyder Equal Area (ISEA) aperture 3 hexagonal cell format [49] with each cell having an area of 32 km 2 and intercell spacing of 6.1 km. A total of 4,235,365 hexagonal cells were used to tessellate the Earth's land surface, which corresponds to an area of 135.5 million km 2 . The global land size excluding Antarctica varies between 134.1 million km 2 and 135 million km 2 based on the scientific literature [50][51][52]. The existing spatial datasets were converted into hexagonal DGGS cells following the approach presented in [53]. The temporal resolution in the research study was determined to be 10 years, and the model was implemented and evaluated using datasets for the years 2000, 2010, and 2020.

Susceptibility Analysis
General multi-criteria evaluation (MCE) approaches [36] have been adopted to implement deforestation susceptibility analysis and were executed at the global level in this research study. Driving factors were identified to represent relevant criteria that characterize the process of deforestation, and susceptibility functions were derived for each criterion. Susceptibility functions transform criterion values into a normalized range between 0 and 1, where 1 indicates highest satisfaction and 0 denotes no satisfaction for the particular criterion. Moreover, each criterion is normalized with the respective suitability function and then weighted and aggregated to obtain the deforestation susceptibility scores for each hexagonal cell. Finally, susceptibility scores can be used to generate global deforestation susceptibility maps that can be one of the inputs that guide the transition rules of the SGA model.
The selected criteria that express some of the key drivers of the deforestation process at the global level are based on the scientific literature and can be grouped into three categories: socioeconomic (population density), terrain (slope, elevation), and proximity (proximity to urban areas, major roads, water bodies, agriculture, forest edge, past forest disturbances) [5,23,54]. Table 1 presents the selected criteria and their respective susceptibility functions as graphs. cells representing forest are converted to the dominant non-forest land-cover type based on the cell's neighbourhood, susceptibility value, and constraint parameter.

Datasets
The study area in this research study encompasses the entire global land surface except for Antarctica, and several spatial datasets with global extent were acquired to implement the model. Land-cover datasets were obtained from the European Space Agency (ESA) portal [43], global roads dataset from the Global Roads Inventory Project (GRIP) [44], protected areas from the World Database on Protected Areas (WDPA) [45], past forest disturbance from the Global Wildfire Information System (GWIS) [46], population density from the LandScan portal [47], and elevation dataset from the United States Geological Survey (USGS) portal [48]. All spatial datasets were converted into Icosahedral Snyder Equal Area (ISEA) aperture 3 hexagonal cell format [49] with each cell having an area of 32 km 2 and intercell spacing of 6.1 km. A total of 4,235,365 hexagonal cells were used to tessellate the Earth's land surface, which corresponds to an area of 135.5 million km 2 . The global land size excluding Antarctica varies between 134.1 million km 2 and 135 million km 2 based on the scientific literature [50][51][52]. The existing spatial datasets were converted into hexagonal DGGS cells following the approach presented in [53]. The temporal resolution in the research study was determined to be 10 years, and the model was implemented and evaluated using datasets for the years 2000, 2010, and 2020.

Susceptibility Analysis
General multi-criteria evaluation (MCE) approaches [36] have been adopted to implement deforestation susceptibility analysis and were executed at the global level in this research study. Driving factors were identified to represent relevant criteria that characterize the process of deforestation, and susceptibility functions were derived for each criterion. Susceptibility functions transform criterion values into a normalized range between 0 and 1, where 1 indicates highest satisfaction and 0 denotes no satisfaction for the particular criterion. Moreover, each criterion is normalized with the respective suitability function and then weighted and aggregated to obtain the deforestation susceptibility scores for each hexagonal cell. Finally, susceptibility scores can be used to generate global deforestation susceptibility maps that can be one of the inputs that guide the transition rules of the SGA model.
The selected criteria that express some of the key drivers of the deforestation process at the global level are based on the scientific literature and can be grouped into three categories: socioeconomic (population density), terrain (slope, elevation), and proximity (proximity to urban areas, major roads, water bodies, agriculture, forest edge, past forest disturbances) [5,23,54]. Table 1 presents the selected criteria and their respective susceptibility functions as graphs. Table 1. Selected criteria of deforestation with their respective susceptibility functions with rationale and criteria weights.

Category
Deforestation Criteria

Susceptibility Functions Rationale Criteria Weights
Socioeconomic Population density Population density is an indicator for the concentration of human activities and thus closeness of possible deforestation processes.

0.03418
Population density is an indicator for the concentration of human activities and thus closeness of possible deforestation processes. Forest areas closer to agricultural land are more prone to deforestation due to expansion of farmlands and agroforestry.

0.12282
Forest areas closer to agricultural land are more prone to deforestation due to expansion of farmlands and agroforestry. 0.12282

0.17331
The susceptibility functions were generated for each criterion and informed by the literature [23,[55][56][57]. The socioeconomic group of criteria rooted in anthropogenic activities is a major determinant of deforestation, and population density is often used as an indicator for the concentration of urban regions and thus human activities [58]. Increasing population density and urban area expansion cause pressure on nearby forests due to the harvesting of wood for construction and fuel, farming, cattle grazing, and urban and infrastructure development [59]. The population density susceptibility function is expressed as a linear membership based on the maximum population density value obtained from the datasets which is 1,168,691 inhabitants per cell. Characteristics of the terrain build another group of criteria. Differences in elevation and slope can represent restrictions to deforestation. Areas with steep slopes are less prone to deforestation as they are unfavourable for other land-use types such as agriculture, infrastructure, and urban development [54,60]. Flat areas, however, allow for accessibility for clearing of forest for agricultural activities and urban and infrastructure development. For the slope susceptibility function, gradients less than 5° yield maximum satisfaction, and slopes steeper than 59° are considered unsuitable for deforestation. The susceptibility function for elevation has maximum satisfaction in locations where altitude is less than 200 m and decreases with increasing elevation until 4900 m. The maximum elevation is set to 4900 m given the highest forest stand is located at this altitude [61,62].
Proximity-based criteria were selected due to different land-use and land-cover features that are driving the deforestation process. Urbanization creates increased demand for land, and deforestation is more likely in areas closer to urban centres as forests are more likely to be cleared for urban expansion, wood harvesting for fuel, and agricultural activities. The function uses a linear membership with maximum susceptibility in locations within 6.1 km of urban areas, which is equivalent to 1 hexagonal cell, and no susceptibility beyond 61 km of urban areas, corresponding to 10 hexagonal cells in the spatial dataset. Major roads provide accessibility to areas dominated by forest land-cover for anthropogenic activities such as urbanization, infrastructure development, agriculture, and resource extraction [63]. Several research studies have indicated a strong positive correlation between deforestation rates and proximity to major roads [64,65]. Prior studies indicate 95% of deforestation can occur with 4.5 km of major roads, with the influence of roads on deforestation extending as far as 100 km [63]. The susceptibility function decreases with increasing distance from major roads, with no susceptibility beyond 97.6 km of roads, which is equivalent to 16 hexagonal cells in the spatial data layer. Proximity to water bodies can also potentially influence the dynamics of deforestation by increasing accessibility to remote areas, transporting forest resources such as timber, and providing

0.17331
The susceptibility functions were generated for each criterion and informed by the literature [23,[55][56][57]. The socioeconomic group of criteria rooted in anthropogenic activities is a major determinant of deforestation, and population density is often used as an indicator for the concentration of urban regions and thus human activities [58]. Increasing population density and urban area expansion cause pressure on nearby forests due to the harvesting of wood for construction and fuel, farming, cattle grazing, and urban and infrastructure development [59]. The population density susceptibility function is expressed as a linear membership based on the maximum population density value obtained from the datasets which is 1,168,691 inhabitants per cell. Characteristics of the terrain build another group of criteria. Differences in elevation and slope can represent restrictions to deforestation. Areas with steep slopes are less prone to deforestation as they are unfavourable for other land-use types such as agriculture, infrastructure, and urban development [54,60]. Flat areas, however, allow for accessibility for clearing of forest for agricultural activities and urban and infrastructure development. For the slope susceptibility function, gradients less than 5° yield maximum satisfaction, and slopes steeper than 59° are considered unsuitable for deforestation. The susceptibility function for elevation has maximum satisfaction in locations where altitude is less than 200 m and decreases with increasing elevation until 4900 m. The maximum elevation is set to 4900 m given the highest forest stand is located at this altitude [61,62].
Proximity-based criteria were selected due to different land-use and land-cover features that are driving the deforestation process. Urbanization creates increased demand for land, and deforestation is more likely in areas closer to urban centres as forests are more likely to be cleared for urban expansion, wood harvesting for fuel, and agricultural activities. The function uses a linear membership with maximum susceptibility in locations within 6.1 km of urban areas, which is equivalent to 1 hexagonal cell, and no susceptibility beyond 61 km of urban areas, corresponding to 10 hexagonal cells in the spatial dataset. Major roads provide accessibility to areas dominated by forest land-cover for anthropogenic activities such as urbanization, infrastructure development, agriculture, and resource extraction [63]. Several research studies have indicated a strong positive correlation between deforestation rates and proximity to major roads [64,65]. Prior studies indicate 95% of deforestation can occur with 4.5 km of major roads, with the influence of roads on deforestation extending as far as 100 km [63]. The susceptibility function decreases with increasing distance from major roads, with no susceptibility beyond 97.6 km of roads, which is equivalent to 16 hexagonal cells in the spatial data layer. Proximity to water bodies can also potentially influence the dynamics of deforestation by increasing accessibility to remote areas, transporting forest resources such as timber, and providing Past forest disturbance areas are often precursors to future forest degradation or agricultural use.

0.17331
The susceptibility functions were generated for each criterion and informed by the literature [23,[55][56][57]. The socioeconomic group of criteria rooted in anthropogenic activities is a major determinant of deforestation, and population density is often used as an indicator for the concentration of urban regions and thus human activities [58]. Increasing population density and urban area expansion cause pressure on nearby forests due to the harvesting of wood for construction and fuel, farming, cattle grazing, and urban and infrastructure development [59]. The population density susceptibility function is expressed as a linear membership based on the maximum population density value obtained from the datasets which is 1,168,691 inhabitants per cell. Characteristics of the terrain build another group of criteria. Differences in elevation and slope can represent restrictions to deforestation. Areas with steep slopes are less prone to deforestation as they are unfavourable for other land-use types such as agriculture, infrastructure, and urban development [54,60]. Flat areas, however, allow for accessibility for clearing of forest for agricultural activities and urban and infrastructure development. For the slope susceptibility function, gradients less than 5 • yield maximum satisfaction, and slopes steeper than 59 • are considered unsuitable for deforestation. The susceptibility function for elevation has maximum satisfaction in locations where altitude is less than 200 m and decreases with increasing elevation until 4900 m. The maximum elevation is set to 4900 m given the highest forest stand is located at this altitude [61,62].
Proximity-based criteria were selected due to different land-use and land-cover features that are driving the deforestation process. Urbanization creates increased demand for land, and deforestation is more likely in areas closer to urban centres as forests are more likely to be cleared for urban expansion, wood harvesting for fuel, and agricultural activities. The function uses a linear membership with maximum susceptibility in locations within 6.1 km of urban areas, which is equivalent to 1 hexagonal cell, and no susceptibility beyond 61 km of urban areas, corresponding to 10 hexagonal cells in the spatial dataset. Major roads provide accessibility to areas dominated by forest land-cover for anthropogenic activities such as urbanization, infrastructure development, agriculture, and resource extraction [63]. Several research studies have indicated a strong positive correlation between deforestation rates and proximity to major roads [64,65]. Prior studies indicate 95% of deforestation can occur with 4.5 km of major roads, with the influence of roads on deforestation extending as far as 100 km [63]. The susceptibility function decreases with increasing distance from major roads, with no susceptibility beyond 97.6 km of roads, which is equivalent to 16 hexagonal cells in the spatial data layer. Proximity to water bodies can also potentially influence the dynamics of deforestation by increasing accessibility to remote areas, transporting forest resources such as timber, and providing water resources for human settlements [66]. The susceptibility function of proximity to water bodies also decreases with increasing distance from water bodies, with no susceptibility beyond 42.7 km of water bodies, which is represented by 7 hexagonal rings. Further, prior studies have indicated agricultural expansion to be another significant driving factor of deforestation [67][68][69]. Forest areas closer to agricultural lands are more likely to be converted for agricultural purposes to support increasing demand for food, biofuels, and animal production. The susceptibility function is based on a decreasing linear function with no susceptibility past 30.5 km, corresponding to 5 hexagonal cells.
Deforestation typically proceeds from the edge of forests into the interior and subsequently leads to the fragmentation of large forest regions into smaller non-contiguous areas [70,71]. Thus, areas closers to forest edges are more prone to deforestation due to their accessibility by the local population. The susceptibility function for proximity to the forest edge decreases with distance until 61 km. Also, past forest disturbance is often seen as a precursor to future forest degradation, and several studies have indicated deforestation is more likely to occur in areas that have experienced some form of disturbance such as forest fire, logging, or mining [72][73][74]. The past disturbance susceptibility function also uses a linear membership where susceptibility decreases with distance until 67.1 km, which is equal to 11 hexagonal cells.

Criterion Weight Generation and Global Susceptibility Maps
Weights were generated for each criterion to reflect the relative importance of the selected criteria in determining deforestation. While criteria weights in most GIS-based MCE methods can be determined by subject experts or stakeholders, this was not possible in this research study due to the global scope of the model's application and resource limitations. Thus, the Automatic Weight Selection technique [75] was applied to generate the weights of importance, and values were normalized so the sum of all weights is equal to one. The technique is based on the comparison between locations of deforestation and random sampling sites using the Cohen's d metric [76]. The obtained criterion weights are presented in Table 1. Based on the normalized criterion values and criterion weights, the Weighted Linear Combination (WLC) technique [77] was used to calculate the overall deforestation susceptibility scores for each hexagonal cell for deforestation. The obtained susceptibility scores were classified into five classes using the equal interval method. They were categorized as very low (0-0.2), low (0.21-0.4), medium (0.41-0.6), high (0.61-0.8), and very high (0.81-1) susceptibility to deforestation. The equal interval method was chosen to classify the susceptibility values due to its ability to create categories of equal sizes and for easy comparison. Figure 2 depicts the obtained global deforestation susceptibility output maps for different parts of the Earth.

Criterion Weight Generation and Global Susceptibility Maps
Weights were generated for each criterion to reflect the relative importance of the selected criteria in determining deforestation. While criteria weights in most GIS-based MCE methods can be determined by subject experts or stakeholders, this was not possible in this research study due to the global scope of the model's application and resource limitations. Thus, the Automatic Weight Selection technique [75] was applied to generate the weights of importance, and values were normalized so the sum of all weights is equal to one. The technique is based on the comparison between locations of deforestation and random sampling sites using the Cohen's d metric [76]. The obtained criterion weights are presented in Table 1. Based on the normalized criterion values and criterion weights, the Weighted Linear Combination (WLC) technique [77] was used to calculate the overall deforestation susceptibility scores for each hexagonal cell for deforestation. The obtained susceptibility scores were classified into five classes using the equal interval method. They were categorized as very low (0-0.2), low (0.21-0.4), medium (0.41-0.6), high (0.61-0.8), and very high (0.81-1) susceptibility to deforestation. The equal interval method was chosen to classify the susceptibility values due to its ability to create categories of equal sizes and for easy comparison. Figure 2 depicts the obtained global deforestation susceptibility output maps for different parts of the Earth.

Deforestation Scenarios
In this research study, three deforestation scenarios were designed and implemented to simulate forest land-cover change under different conditions: Business as Usual (BAU), Accelerated Deforestation (AD), and Sustainable Deforestation (SD) scenarios. The Business as Usual scenario assumes the historical rate of deforestation observed between 2010 and 2020 will continue in the future [78]. Further, deforestation inside protected areas is allowed due to ineffective implementation of forest conservation policies in some regions and to allow either for urbanization or agricultural expansion [79]. Research reveals that increased demand for forest products and services could amplify rates of deforestation by half in some parts of the world [80]. Therefore, the Accelerated Deforestation scenario assumes the rate of forest loss would be 50% higher than the current rate as well as loose implementation of environmental conservation policies. This represents a pessimistic deforestation scenario used to characterize forest-cover change under complete absence of forest management policies and lack of political commitments to reducing deforestation at the local and global levels. The Sustainable Deforestation scenario is the most optimistic and assumes a reduction in the current trend of deforestation by 50% through 2050 and 75% by 2100. Also, it encompasses the strict enforcement of forest conservation policies and programmes, and deforestation is restricted in protected areas under this scenario. Prior studies indicate the effective implementation of forest conservation policies and measures can positively impact global climate change. The research findings indicate reducing rates of deforestation by 50% could potentially reduce carbon emissions from land-cover change by 13 to 50 gigatons of carbon (GtC) [81,82].
Under each scenario, the model is constrained by the rate of deforestation at the country level and calculated using the 2010 and 2020 land-cover datasets. However, due to a lack of adequate data, the conversion of forest cover to water and snow and the process of reforestation are not considered.

Model Implementation and Evaluation
The spherical deforestation model was implemented in the Python programming language [83] using the DDGRID open-source library [84]. The model was implemented on a workstation with an Intel(R) Xeon(R) Gold 6128 CPU @ 3.40 GHz 3.39 GHz processor and 32 GB RAM, with the processing time for each scenario implementation varying between 81 and 88 h. The model was run for eight iterations with a temporal resolution of 10 years to simulate global deforestation between 2020 T i and 2100 (T i+8 ) and for the three scenarios.
Model evaluation was performed using the relative operating characteristic (ROC) technique [85] and the Figure of Merit (FoM) [86]. ROC entails metrics for assessing the performance of binary classification with continuous output or rank order values [87,88]. ROC applies thresholds to generate a contingency table with four performance descriptors: true positives (TP), false negatives (FN), false positives (FP), and true negatives (TN) [89]. True positives (TP) correspond to changed forest cells correctly simulated as change by the model, false negatives (FN) are unchanged forest cells wrongly simulated as changed cells, false positives (FP) are changed forest cells the model was unable to simulate as changed cells, and true negatives (TN) are changed forest cells simulated as change but to the wrong land-cover class [90]. From the contingency table, the true positive rate (TPR) and false positive rate (FPR) can be calculated as follows: By plotting the TPR on the vertical axis and FPR on the horizontal axis of the graph, the ROC curve and Area Under the Curve (AUC) metrics are obtained [91]. AUC values range between 0 and 1, where a larger value indicates higher model accuracy. Additionally, the simulation outputs are compared with actual land-cover datasets using the Figure of Merit (FoM) index which can be expressed as follows: where the definition of the term misses is the same as false negatives, hits denote true positives, and false alarms are equivalent to false positives. In this research study, global datasets for the period 2000-2010 were used for model calibration, and then global datasets for the period 2010-2020 were used for model validation. The AUC value obtained for the global spherical geographic deforestation model was 0.9 in the calibration phase and 0.87 in the validation phase. Other global geosimulation model applications in the scientific literature report AUC values ranging between 0.72 and 0.93 [92][93][94]. For the FoM metric, the value obtained was 36.5 for the model calibration and 29.9% during model validation. FoM values obtained in other global geosimulation models range between 19% and 43% [33,93]. Thus, the evaluation of the proposed modelling methodology yields commensurate values.

Global and Regional Variations in Forest-Cover Change
In 2020, forest covered 48 million km 2 of the terrestrial Earth's surface, corresponding to 35.6% of the global land area. The simulation results of deforestation are presented in Figure 3 for North America only as an illustration of detailed model outputs, and for each time step between 2020 and 2100 under the Accelerated Deforestation (AD) scenario. The obtained simulation outputs of deforestation under the different scenarios by the year 2100 compared with the base year 2020 are also presented for different parts of the globe in Figure 4.
Under the BAU scenario, the global forest extent shrinks to 43 million km 2 , decreasing by 5 million km 2 between 2020 and 2100, which corresponds to an annual forest loss of 63 thousand km 2 per year, about twice the size of the Netherlands. Approximately 10.5% of the forest extent in 2020 would be lost by the end of the 21st century based on the current trend of deforestation. Under the AD scenario, the global forest extent would decrease by 7.3 million km 2 , representing a forest loss of 15.3% by 2100. Conversely, the simulation results indicate 1.8 million km 2 of the global forest area would be deforested by 2100 under the SD scenario, which represents a global forest loss of 3.7%. Figure 5 presents the global cumulative deforestation obtained for the three scenarios.      The simulation results also revealed marked differences in deforestation dynamics at the continental level. Table 2 presents summaries of the simulated forest-cover change per continent in the period between 2020 and 2100. At the continental level, Europe had the largest forest-cover loss with 1.9 million km 2 , corresponding to 15.5% of the continent's forest extent in 2020, under the BAU scenario. The results further indicate forest loss in Europe would reach 2.7 million km 2 under the AD scenario. It must however be noted that the Russian Federation which is considered part of Europe in this research accounts for 81.8% of Europe's forest area in 2020. The simulation outputs also reveal considerable deforestation in North America with 1.02 million km 2 of forest lost under the BAU scenario and 1.5 million km 2 under the AD scenario. Forest losses in Africa, Asia, and South America are revealed to be 0.58 million km 2 , 0.66 million km 2 , and 0.76 million km 2 , respectively, under the BAU scenario. Figure 6 depicts simulated deforestation in 2100 compared to the base year for different forest regions across the globe under the three scenarios. The Amazon Forest, Congo Basin, and Eastern USA regions are presented here due to the extensive forest loss revealed by the simulation results. Considerable differences in forest extent and rates of deforestation can be observed at the country level as well. The simulation results indicate the largest extent of deforestation would occur in the Russian Federation, Canada, the United States, and Brazil. Under the AD scenario, for example, 2.3 million km 2 , 0.81 million km 2 , 0.67 million km 2 , and 0.45 million km 2 of forest area were lost in the Russian Federation, Canada, the United States, and Brazil, respectively, between 2020 and 2100. In relative terms, 17.4% of the forest in Canada and 20.3% of the forest in the United States would be lost by the end of the 21st century under this scenario. Moreover, there are several countries (over 50) with zero forest loss due to these countries having no forest cover, a lack of data, or the country being too small to capture the change in forest-cover change. Considerable differences in forest extent and rates of deforestation can be observed at the country level as well. The simulation results indicate the largest extent of deforestation would occur in the Russian Federation, Canada, the United States, and Brazil. Under

Forest Change in Protected Areas
In 2020, about 9.3 million km 2 of the global forest was in protected areas, corresponding to 19.5% of the global forest extent. Simulation results indicate that this proportion, however, increases to 21.8%, 23%, and 20.2% by 2100 under the BAU, AD, and SD scenarios, respectively. The rate of deforestation within protected areas was significantly lower than the global rate as well. Between 2020 and 2100, only 0.03% of forests in protected areas were deforested under the BAU scenario, and 0.09% of forests in protected areas were deforested under the AD scenario. Conversely, the forest extent in protected areas under the SD scenario remained the same over the simulation period with no forest loss due to the restriction of deforestation in these areas. Table 3 presents the percentages of the share of forest extent located in protected areas at the continental level by 2100 and compared to the base year under the three scenarios. By 2100, 3.7 million km 2 of forest in South America would be in protected areas under the BAU scenario, which represents 38.7% of all forest cover on the continent. This represents the largest share of protected forest area at the continental scale. In contrast, Oceania, North America, and Europe have the lowest percentage of protected forests at the continental level under all scenarios.

Discussion
Based on the simulation outputs, the global forest extent is projected to decrease over the coming decades, with the total forest extent in 2100 ranging between 46 and 40 million km 2 . However, the rate and magnitude of deforestation differ among the three scenarios and at the regional and country levels. The results indicate that Europe has the largest extent of deforestation by 2100; however, 85.4% of the deforested areas in Europe would occur in the Russian Federation. The results from the model presented in this research are comparable with other projections found in the literature with a range of outcomes varying between 25 million and 50 million km 2 by 2100 under different scenarios [30,31,95]. However, as a caveat, simulation results from the presented model are largely dependent on the quality of the datasets utilized in the research. For instance, the rates of deforestation obtained at the country level for implementing the scenarios are derived from the available ESA-CCI land-cover datasets.
The pattern of forest-cover change observed also follows the process of deforestation as reported by prior studies [70,71]. From the simulation outputs, deforestation initially begins from the fringes of large forest regions and diminishes into the interior. This spatial pattern can be justified where forest regions in proximity to past forest disturbances, urban areas, water bodies, and road networks are more prone to deforestation. Over the course of the simulation run, regions initially covered by large forests become fragmented as deforestation spreads into the forest core. This pattern can be observed from the temporal simulation outputs presented in Figure 3.
While no forest management policies were explicitly included in the scenario design and implementation, the difference in results among the scenarios indicates the model can be utilized to assess different forest conservation policies. When properly implemented, protected areas can be used as effective forest conservation schemes to reduce deforestation. For instance, in the SD scenario where forest management policies are strictly implemented and deforestation in protected areas is not allowed, the simulation results reveal no loss of forest cover. However, such an outcome would be difficult to achieve due to the financial and human resources required to implement such a policy across large regions. However, the positive impact of protected areas on deforestation is proven by other studies in the scientific literature [63,96]. A research study [97] indicates that deforestation in protected areas accounted for only 5% of net forest loss between 2000 and 2021 in the Brazilian Amazon. Presently, protected areas are predominantly located in tropical regions of Africa, Asia, and South America, in contrast to the small proportion of protected forest areas in Europe and North America as obtained from the simulation results. According to [98], less than 10% of subtropical humid forests, temperate steppe and boreal coniferous forests largely found in Europe and North America are protected.
Despite the model's capabilities in simulating deforestation at the global level and across different regions, the research study has some limitations. Improvements can be made to the model depending on the availability of quality and detailed global datasets to simulate land-cover change at finer spatial resolution as well as including more relevant criteria related to deforestation such as soil properties, climate variables, and forest type classification. By experimenting with different weighting methods such as Analytic Hierarchy Process (AHP) [99] or incorporating advanced spatial decision techniques such as Ordered Weighted Averaging (OWA) [100] and Logic Scoring of Preference (LSP) [101], the deforestation susceptibility analysis can be further improved for the SGA modelling framework. The selection of the relevant drivers of deforestation and generation of criterion weights can also be determined through engagement with subject experts and stakeholders in the model implementation phase of the research study. Considering the spatial heterogeneity and dynamics of the different drivers of deforestation across diverse regions and countries, implementation of region-specific economic or climate policies would be beneficial to improve the model. While only homogeneous forest was considered in this research, different forest types such as rainforest, boreal, deciduous, and mangrove, to name a few, can also be included in order to incorporate detailed characteristics of the different forest change dynamics. Additionally, the model can further be developed to incorporate multiple land-use/land-cover change types to reflect their different dynamics. This can assist in making informed decisions at country or regional levels and considering REDD+ [102] and OECD [103] concepts. Consideration of the natural regeneration of forests, reforestation, afforestation, and age of forested areas would be beneficial to enhance the proposed modelling approach. Moreover, with climate being one of the drivers of forest distribution, the inclusion of climate variables and scenarios can enhance the model's ability to characterize future deforestation patterns and consider the effects of climate change. Augmenting computational power with more efficient code for the SGA model to run faster or on multiple processors when using global datasets would be another advantage.

Conclusions
This research study aimed to develop and implement a unique spherical geographic automata modelling approach that has been applied to represent and simulate the global deforestation process. Compared to existing geosimulation models and applications, the proposed model considers the curvature of the Earth's surface, which is often ignored when modelling at the global level. For large-scale spatial applications, it is determined that the use of planar spatial models can produce very different results. Further, most global land-cover change models in the literature generally focus on simulating agricultural and urban land-use change, with few studies including deforestation.
Results from this research study indicate the spherical deforestation model can be successfully implemented to simulate the forest land-cover change process at the global level and under different "what-if" scenarios. Ultimately, the model is flexible and allows for further enhancements. Thus, the proposed deforestation modelling approach has a solid foundation to be used by intergovernmental entities, policymakers, ecologists, and researchers to support global forest management and conservation. With the United Nations (UN) Sustainable Development Goal (SDG) 15 set on promoting the sustainable use of terrestrial ecosystems and sustainable forest management and halting biodiversity loss, the spherical geographic deforestation model proposed in this research study can provide valuable insight and a spatial decision-support tool for achieving these targets.