Geospatial Simulation Model of Deforestation and Reforestation Using Multicriteria Evaluation

: Deforestation is an anthropic phenomenon that negatively a ﬀ ects the environment and therefore the climate, the carbon cycle, biodiversity and the sustainability of agriculture and drinking water sources. Deforestation is counteracted by reforestation processes, which is caused by the natural regeneration of forests or by the establishment of plantations. The present research is focused on generating a simulation model to predict the deforestation and reforestation for 2030 and 2050 using geospatial analysis techniques and multicriteria evaluation. The case study is the North Paciﬁc Basin, which is one of the areas with the greatest loss of forest cover in Mexico. The results of the spatial analysis of forest dynamics determined that the forest area in 2030 would be 98,713.52 km 2 , while in 2050 would be 101,239.8 km 2 . The mean annual deforestation and reforestation expected in the study area is 115 and 193.84 km 2 , for the 2014–2030 period, while mean annual deforestation and reforestation values of 95 and 221.31 km 2 are expected for the 2030–2050 period. Therefore, considering the forest cover predicted by the deforestation and reforestation model, a carbon capture of 16,209.67 ton / C was estimated for the 2014–2030 period and 587,596.01 ton / C for the 2030–2050.


Introduction
Over the years, the land surface has undergone important transformations that have led to climate change, with deforestation being one of the most prevalent problems in the loss of biodiversity [1][2][3][4], alteration of hydrological cycles [5] and the loss of carbon (C) stocks [6]. Deforestation causes deterioration of ecosystems and the services they provide to human beings, such as drinking water, medicines, food, fuels and, of course, climate regulation [7][8][9].
Given this problem, strategies must be developed to counteract deforestation processes effectively since forest cover is very important. Forests must be considered biological systems with multiple purposes, such as storing the greatest amount of carbon and reducing carbon dioxide (CO 2 ) emissions [10,11]. Therefore, the management of forest resources must consider different aspects related to the conservation and sustainability of the ecosystem to avoid degradation and to preserve biological species.
For this reason, the Reduction in Emissions from Deforestation and Forest Degradation (REDD+) strategies have been implemented, which require spatial information about the changes in forest cover and the development of models to project the forest processes in the future, and consequently the CO 2 emissions associated with the greenhouse effect. For that purpose, the use of geospatial models is considered as an alternative due to the fact that they provide knowledge on the dynamic processes in a geographic space [12]. By using the geospatial models, alternative scenarios and strategies can be developed, where the impact and potential consequences of the deforestation and reforestation strategies can be evaluated [13,14].
Several geosimulation studies have been proposed to describe forestry processes using different techniques, such as cellular automata [15][16][17][18]; multiagent systems [19][20][21]; integrated multiagent systems and cellular automata [22]; Markov chain [23]; Multicriteria Evaluation (MCE) techniques [24,25]; Bayesian networks [26]; conversion models of land use and its effects [27], as well as a platform called ForestSim based on Agent-Based Modeling was designed that integrates tools to model forest processes [28] and similarly, a geocomputational methodology has been also proposed based on a geographic information system and decision support system [29]. Most of these models simulated forest processes based on biophysical, socioeconomic and accessibility factors, and were developed to describe the anthropogenic impacts and to measure strategies to mitigate the loss and deterioration processes of forest covers [30].
In Mexico, some simulation scenarios have been proposed for forests. One of the most representative studies is a map of deforestation risks, which was generated by using linear regression between the deforestation rate and some spatial variables [31]. Another example is the research developed by [32], in which a simulation model was developed using MCE and Geographic Information Systems (GIS) techniques with the aim of locating suitable areas for forest plantations in the state of Mexico.
Similarly, future scenarios for 2025 were modeled using cellular automata in the Sierra Madre Oriental in San Luis Potosí [33]. In Mas and Flamenco [34], two scenarios were generated using cellular automata and weights of evidence techniques in a region of southeastern Mexico: the trend-oriented and the sustainable scenarios. The land-use cover scenarios showed a decrease in the exchange rate and a relocation of forest clearings in secondary areas. Likewise, future scenarios of changes in forest cover were defined for 2025 in the state of Oaxaca using cellular automata and weights of evidence techniques [35].
Most of the studies that have been carried out in different parts of the world had been focused on determining areas to be deforested [26,[36][37][38] or reforested [39,40], without considering the occurrence of both forestry processes. Recent studies are focused on the genetic programming of forest areas [41], but the impact of forest cover changes on carbon sequestration has not been exhaustively analyzed. Some studies have presented numerical results of the carbon loss and capture by deforestation and reforestation processes [42][43][44]. However, a spatial analysis of land use and cover related to carbon loss and capture has been not discussed. The carbon capture and loss that could be produced by forest processes in the future has not been addressed either.
The main objectives of the present study are: (1) to analyze the forest dynamics in the North Pacific basin, Mexico, using geospatial analysis; (2) to generate a geospatial simulation model by using multicriteria evaluation techniques that includes the future behavior of the deforestation and reforestation for 2030 and 2050. Since forest territorial images were generated to describe deforestation and reforestation processes in the study area, the present study cannot differentiate between natural forest regeneration and forest plantations. Based on these results, the carbon capture and loss were estimated and mapped.

Study Area
The North Pacific basin is located within the entire state of Sinaloa and also comprises some municipalities of the states of Durango, Chihuahua, Zacatecas and Nayarit in Mexico. The basin has an area of 152,013 km 2 , of which 64.10% corresponds to forests (primary and secondary coniferous, oak, mountain mesophilic, deciduous, evergreen, sub-deciduous forest, and cultivated forest) ( Figure 1). The North Pacific basin represents 8.0% of the surface of Mexico according to the National Institute of Statistics and Geography (INEGI) [45]. The estimated population in the basin is 4,466,000 inhabitants [46]. The population increase and their socioeconomic activities have caused changes in the occupation and land use in some regions of the basin, which have produced the loss of forest cover and therefore negative effects on ecosystems, such as soil degradation, decreases in aquifer storage due to changes in the water cycle and loss of biodiversity, as stated by the National Water Commission (CONAGUA) [47]. As previously described, deforestation has a direct effect in the study area because it is the most fertile valley in Mexico, contributing 30.2% of the cultivated land [48]. Agricultural production is mainly due to the presence of many surface water bodies in this study area, where 13 major rivers and 11 large dams are located throughout the basin and serve different functions, such as water storage, fishing, and irrigation [47].

Data
Several variables, for example altitude, slope, forest production, restoration areas, soil types, population density, and degraded forest, among others, were considered for the geospatial simulation model ( Table 1). The selection of these variables was based on an expert consultation to describe the possible causes of deforestation [49] and was complemented with an extensive bibliographic review [15,31,50,51] and an exhaustive assessment that took into account the availability of information, the objectives of the study and the scenarios to be simulated. Likewise, the official Land-use and vegetation maps of the INEGI at a scale of 1:250,000 were used to perform the analysis of the forest cover [45], and the National Greenhouse Gas Inventory database of Mexico reported by [52].

Methodology
The geospatial simulation model of deforestation and reforestation was established using the forest cover of 2014 as baseline. The MCE methods were integrated to model deforestation and reforestation of forests. The implementation of MCE techniques began with the definition of the simulation model objectives, as part of the design of the geospatial model. Then, some criteria, factors and restrictions were established for the selection and preprocessing of the data. Finally, the land use demand was calculated and simulation maps for deforestation and reforestation were depicted. The methods and techniques that were used to simulate the geospatial model are described in greater detail in the following section.

Selection and Preprocessing of Data
First, the official maps of Land-Use and Vegetation (LUV) of the INEGI [45] were used to obtain the analysis of deforestation and reforestation. A geometric and topological correction was applied, and the categories of LUV maps were homogenized according to the Biennial Update Reports (BURs) of the National Greenhouse Gas Inventory of Mexico reported to the United Nations Framework Convention on Climate Change (UNFCCC) [52].
The variables selected for the simulation of the geospatial model were: population density, marginalization index, distance to localities with less than 2500 inhabitants, soil types, distance to roads, distance to hydrographic features, distance to agriculture, distance to grasslands, distance to protected natural areas, forest production, altitude, slope, areas for forest restoration and degraded forests. The variables were spatialized (alphanumeric variables), standardized (legend of classes and uses) and rasterized (vector variables) with a pixel size of 100 m.

Demand for Land Use
A land use demand model was developed to simulate the forest area behavior in the future. The Land-use demand simulation consisted of estimating the forest area that would be covered in 2030 and 2050 based on a temporal analysis of the forest areas in 1986,1993,2002,2007,2011 and 2014 using a third-order polynomial function, of which the precision and uncertainty were analyzed.
Then, the deforestation and reforestation covers for the future scenarios were defined. The annual deforestation and reforestation rates for the future scenarios were established based on the behavior observed for both forest processes in the periods of 1986-1993, 1993-2002, 2002-2007, 2007-2011, and 2011-2014. The deforestation and reforestation scenarios were in accordance with the goals of the National Forest Program of Mexico [53], the Sustainable Forest Development Plan of Mexico [54], and the International Sustainable Development Objectives (SDO) of the Intergovernmental Panel on Climate Change (IPCC) regarding the current and future global state of forests [55], which are aimed to reduce deforestation rates to more than 50% and increase reforestation for 2030 and 2050.

Geospatial Simulation Using MCE
The Weighted Linear Summation (WLS) was used to obtain the most suitable areas for deforestation and reforestation. The WLS method is based on both weighted factors and constraints [56,57]. Three criteria (biophysical, socioeconomic and proximity) were used. Each criterion was weighted based on external expert opinions and a bibliographic review. The values obtained were later normalized in such a way that the sum of the weights of the three criteria was equal to 1.
Once the criteria were defined and weighted, the variables were modelled according to the objectives of each scenario using land mapping. The biophysical factors (soil types, land use and vegetation, areas for restoration and degraded forests) were characterized by thematic classes. The Saaty analytical hierarchy method was applied with the purpose of establishing a relative hierarchy of importance between the categories of the same factor [58].
The assignment of weights to each thematic class was made considering quantitative and qualitative aspects. From the quantitative point of view, a spatial intersection was made between the variables of deforestation and reforestation during the different periods and the percentage of impact of each class was determined. From the qualitative point of view, the opinions of 8 external experts were used: 3 from the National Forestry Commission (CONAFOR), 3 from the Ministry of the Environment and Natural Resources (SEMARNAT) and 2 from the National Commission of Protected Natural Areas (CONANP). The experts established the weights for the categories of the same factor.
A byte scale of 0-255 was used for factor normalization by using the multi-criteria evaluation module in the TerrSet software. The highest value of the scale represents the highest probability of belonging to the decision group and fuzzy logic techniques were applied according to [59,60]. Table 2 shows the maximum and minimum values observed for each factor. Table 2 also shows the functions and the scale used for the normalization of each factor.
The model restrictions were binary alternatives, where the value of 0 indicated the areas excluded from the analysis and the value of 1 corresponded to the possible locations for simulation. The deforestation model restrictions included all categories of land-use cover, except forest areas where future deforestation is estimated to take place. Regarding the reforestation model, the covers were enabled based on the ability of the areas to be converted into forest areas, such as scrub, pasturelands, areas for restoration, degraded forests, other types of vegetation and other lands. Once each of the weighted and normalized factors was obtained, the WLS equation (Equation (1)) was used to evaluate the different alternatives. The WLS equation consisted of adding the result of multiplying the individual pixel value of each normalized factor by the final weight assigned to that same factor using the following equation: where r i is the adequacy level of alternative i, w i is the weight of factor j, and v ij is the weighted value of alternative i of factor j.
Since future conditions can differ from the historical data, an uncertainty analysis was carried out to describe the range of potential forest cover in the North Pacific Basin at the 95% confidence level. Besides, the Extended Fourier Amplitude Sensitivity Test (E-FAST) technique suggested by [61] was used to determine the stability/robustness of the model results through the first order sensitivity indices (S i ). The E-FAST technique is based on variance estimation and uses a Monte Carlo simulation for generating random samples (probability distributions) for each of the independent input factors used in the geospatial model [62]. The E-FAST test did not require any particular assumptions on the model structure, and provided quantitative sensitivity measures for each factor. The first order sensitivity indices (S i ) measured how the i th factor contributes to the variation of the model results without taking into account the interactions with other input parameters. A detailed description of the E-FAST method is given in [61].
The sensitivity analysis using the E-FAST technique was implemented according to the following procedure: first, the frequency distribution of the factors was established; then, a sample of the different factors used in the model was extracted. The sensitivity analysis consisted in fluctuating each of the input factors. Each of the model factors' values were decreased and increased by 25%, while the others were steady. The model was executed 5100 times. The forest cover deviations were used to calculate sensitivity indices and reflected the model sensitivity. Finally, the sensitivity indices (S i ) were obtained and analyzed.
Once the uncertainty and robustness of the geospatial simulation model was evaluated, the suitability maps (deforestation and reforestation maps) were obtained from the geospatial simulation model. Then, a multiobjective assignment (MOLA: multiobjective land allocation) was carried out with the assumption that the deforestation and reforestation covers are conflicted and interacted with each other [63].

Indicators and Mapping of Carbon Loss and Capture
An intersection between the land-use and vegetation (LUV) map of the baseline year (2014) and the map obtained from the geospatial simulation model was used to determine the indicators of change and sustainability for the simulated deforestation and reforestation processes. The carbon balance was estimated based on the carbon loss by deforestation and carbon capture by reforestation. The carbon balance is related to live aerial biomass and roots and is expressed in tons of Carbon (ton/C), as indicated in the National Greenhouse Gas Inventory in Mexico reported by [52]. Therefore, the distribution maps of carbon loss and capture were obtained for the 2014-2030 and 2014-2050 periods.

Land Use Demand
The analysis of deforestation and reforestation in the study area is presented for during the period of 1986-1993, 1993-2002, 2002-2007, 2007-  Two different ways of reforestation can be recognized in the study area: natural forest regeneration and forestry plantations. However, no differentiation between both processes can be completely identified since only 100 ha are related to forest plantations, according to official information provided by INEGI [45]. Based on deforestation and reforestation analysis, a demand model for the land area covered by forests was developed (Figure 3). The land demand model obtained is a simplification of the real forest cover behavior observed in the North Pacific Basin. The assumptions, input data and other criteria used for the model development are representative of the real conditions and they can be inaccurate. Therefore, measuring the precision with which the demand model can reproduce observed conditions is very important since the usefulness of the model depends on its accuracy and reliability. In the present study, the forest cover time series was adjusted to a 3rd order polynomial function. The equation of the fitted model is y = 0.1986x 3 − 1189.7x 2 + 2, 375, 098.7x − 1, 580, 264, 028.01 The third order polynomial model observed a reliability of 98.8%, which could be considered as a good fit ( Table 3). The performance of the land demand model in fitting the historical data was also assessed by using the standard error of the estimate (SE) and the mean absolute error (MAE). The standard error of the estimate shows that the standard deviation of the residuals is 279.968. This value was used to construct prediction limits for future observations. The MAE of 173.272 represents the average value of the residuals. Since future conditions can differ from the historical data, an uncertainty analysis was carried out to describe the range of potential forest cover in the North Pacific Basin at the 95% confidence level (Figure 4). The range of uncertainty is related to a frequency distribution that characterizes the variability of deforestation and reforestation in the study area. The variability observed under future scenarios can be related to inadequate information, incorrect assumptions and the natural variability of forest processes (in particular, natural variation related to reforestation).
According to model results, a forest cover of 98,713.52 km 2 was estimated for 2030, while a value of 101,239.8 km 2 was estimated for 2050. The deforestation and reforestation areas and their mean annual values for 2030 and 2050 were determined in Table 4.
The results show that the forest cover estimated for 2050 was similar to that observed in 1986. The forest cover recovery estimated for 2050 can be explained by the annual averages of reforestation expected for the 2014-2030 and 2030-2050 periods, which were greater than the annual average deforestation projected.

Geospatial Simulation Using MCE
The weighting criteria results and the factors assigned by the consulting experts and the bibliographic reviews are summarized in Table 5. The weighting criteria (W c ), the weighting factors (W i ), and the final weights assigned to the factors (W j ) are identified.  Table 6 shows the weighting results by the categorical factors: soil types, land use and vegetation, degraded forests and areas for restoration. The categorical factors results were obtained using the Saaty pairwise comparison method, and their importance was assigned to each class according to the opinion of experts, resulting in a consistency coefficient (c). Through fuzzy logic techniques, 16 normalized factors were obtained with a scale from 0 to 255 ( Figure 5). It should be noted that the Land-use and vegetation factor was normalized twice, depending on the objective to be simulated (deforestation and reforestation). Additionally, all other factors were only normalized once for both simulation objectives. The model restrictions for deforestation and reforestation were depicted, respectively ( Figure 6). For the deforestation objective, only the forest categories were activated, while for the reforestation, the categories of scrub, pasturelands, areas for restoration, degraded forests, other types of vegetation and other lands were activated. Once all normalized factors, restrictions and weights were obtained, the WLS equation was applied to generate suitability maps for deforestation and reforestation (Figure 7). According to the results obtained from E-FAST, the factors with the highest weights were considered as the most sensitive factors in the model (Table 7). Any change in any of the most sensitive factors represents a significant impact on the model. The most important factor was Forest Production for both objectives (deforestation and reforestation) with a S i of 0.236 and 0.416, respectively. For the deforestation objective, other factors with greater importance were the Marginalization Index, Distance to Protected Natural Areas, Slope, and Distance to Agriculture with S i indices of 0.114, 0.102, 0.093 and 0.079, respectively. For the reforestation objective, the most sensitive factors were the Slope and Altitude with S i values of 0.183 and 0.073, respectively, followed by the Distance to Agriculture (S i = 0.059) and Soil Types (S i = 0.046). Afterwards, a classification (ranking) of each of the resulting objectives (suitability maps) was carried out, and a multiobjective allocation was made to solve possible conflicts. Thus, the scenarios proposed for 2030 and 2050 were obtained according to the required land demand for each objective (Figure 8).

Indicators and Mapping of Carbon Loss and Capture
To analyze the future perspectives of the model, an impact analysis of the land use and vegetation in the basin was carried out. According to the forecast results, deforestation by 2030 would occur mainly in areas of primary deciduous forest (33.60%), primary coniferous forest (31.87%), secondary deciduous forest (13.85%) and primary oak forest (10.37%); meanwhile, reforestation in 2030 would take place on grasslands (86.77%) and degraded forests (11.64%) (Figure 9). Regarding the 2050 scenario, the deforestation model indicated that the most impacted categories would be primary deciduous forest, primary coniferous forest, secondary deciduous forest and primary oak forest, with deforested area percentages of 36.35%, 29.35%, 14.00% and 10.73%, respectively, regarding the total deforested area expected. The reforestation in 2050 would be concentrated in the grasslands and degraded forest categories, with 83.85% and 13.67%, respectively ( Figure 8).
Subsequently, the indicators of carbon capture and loss in live aerial biomass and roots were estimated with respect to the deforested and reforested covers of the proposed scenarios. For the 2014-2030 period, 105,428.04 ton/C would be lost due to deforestation, while 121,637.71 ton/C would be captured by reforestation. During the 2030-2050 period, it was estimated that 109,242.02 ton/C could be lost, while 696,838.03 ton/C could be captured ( Table 8). The carbon capture and loss scenarios indicate a forest cover recovery in the study area, with a positive impact from the environmental point of view. The results obtained are promising in terms of climate change mitigation at local and regional level.
Finally, the carbon loss and capture distribution maps were generated, which in turn were referred to the geospatial simulation model of 2030 and 2050. The spatial distribution of the carbon loss by deforestation in the different categories was depicted in the same way as the carbon capture by reforestation (Figure 9).

Discussion
The results showed that the present study contributed significantly to the management of forest ecosystems and to decision-making processes by providing information on the future state of forest cover in North Pacific Basin in Mexico.
First, a novel procedure is proposed to estimate the land use demand for the forest cover, by combining a retrospective analysis (time series analysis) and a non-linear regression (3rd order polynomial function) to describe the forest cover behavior towards the years 2030 and 2050. Based on a historical analysis of deforestation and reforestation in the study area, the average annual rates of both forest processes were estimated. This procedure can be considered novel because other simulations usually only use the land use tendency, such as Markov chains [64] and cellular automata [65].
The proposed model could be considered as a sustainable alternative for all the stakeholders and decision-makers in government since future deforestation is proposed on degraded forest areas, which will not affect forests with high environmental importance. Likewise, the expected reforestation showed a positive impact in terms of mitigating the emission of greenhouse gases in the basin.
Most of studies in the literature have focused on a single objective, such as developing deforestation scenarios [26,[36][37][38]66], and few studies have included the reforestation [39,40]. However, no studies were found that consider both forest processes at the same time to generate the future simulation of the forest cover, since conflicts arise when determining whether an area can be deforested or reforested in the future.
Besides which, the proposed model includes variables, such as forest production, degraded forests, and areas for restoration, that have not been included in other studies to determine the most suitable areas for deforestation or reforestation. The weighting criteria and the weights assigned to categorical variables also differ from other studies that have used MCE [32,67,68]. The weights of categorical variables were assigned from two points of view: the qualitative criteria, which was based on the opinion of experts and the quantitative criteria, which considered spatial computations.
Sometimes, the forest classification information provided by official sources can be ambiguous. For instance, an in-depth classification of the mountain forests in the study area can be found in [69], which differs from the classification provided by INEGI [45]. Therefore, the development of geospatial models must be supported with a highly detailed classification of forests (such as biogeographic or ecologic classification), which could be of interest for other specialists.
Natural regeneration of the forest is a normal induced ecological phenomenon of succession, while forest plantation (especially with non-native species) is an ecologically undesirable practice [70][71][72]. Differentiating the natural reforestation process from the induced planting process is a common problem in geosimulation that could be solved using very high resolution data (spatial, spectral and radiometric) derived from remote sensors [73,74]. Therefore, a more in-depth information from an ecological point of view could be obtained through the remote sensing analysis of high resolution data.
The geospatial model estimated the carbon loss and capture by deforestation and reforestation processes. As a result, distribution maps of the carbon loss and capture for 2030 and 2050 were obtained, generating indicators to support the decision-making process for forest management, mitigation and planning purposes, which are directly related to climate change at the local, regional and global levels. Given that some studies suggest that carbon storage comes directly from forest cover [75,76], the present study can also help to visualize future distribution of ecosystem services.
Finally, the proposed methodological approach for the implementation of the scenarios could be useful for generating other ecological indicators and can be applied to other basins considered as a priority by the institutions involved in forest planning and management.

Conclusions
The geospatial simulation model provides relevant information on a possible future state of the forest cover in North Pacific basin. Based on the scenarios simulated for 2050, the forest cover in the study area would be similar to the observed in in 1986, which would be fulfilling the National Forest Program goals.
The suitability areas of deforestation and reforestation were estimated and mapped. This simulation showed a decrease in the mean annual rates of deforestation, where mean annual rates of 115 km 2 and 95 km 2 are expected for the 2014-2030 and 2030-2050 periods. During the same period, an increase in the mean annual reforestation rates was identified, with mean annual rates of 193.84 km 2 and 221.31 km 2 , respectively. The sensitivity analysis using E-FAST identified the most significant factors in the deforestation and reforestation models. In particular, the Biophysical factors were recognized as the most important parameters of the model.
Based on the results obtained in the present study, the geospatial simulation model showed sustainable simulation scenarios from a forest management point of view. The model uncertainty analysis showed that even in the worst case scenario, the North Pacific basin shows a recovery in forest cover by 2050. The loss and carbon capture were estimated and mapped, considering the forest cover predicted by the deforestation and reforestation model. A carbon capture of 16,209.67 ton/C was estimated for the 2014-2030 period and 587,596.01 ton/C for the 2030-2050 period. The most affected land use categories were the primary deciduous forest, primary coniferous forest and secondary deciduous forest. In contrast, grassland and degraded forest categories were identified as the most suitable covers for reforestation. This information is of a great importance for government institutions for the management and planning of forest ecosystems.
The present study demonstrated that the geospatial analysis techniques and multicriteria evaluation are viable strategies for the development of forest management scenarios, since they take into account different perspectives such as the experts' opinions, national and international goals, and a set of biophysical, socioeconomic and proximity variables that generated a robust methodology that can be used to predict deforestation and reforestation processes in different regions around the world.