Cuadernos de Investigación Geográfica / Gegographical Research Letters no 44 (2, 2018)

The main objective of this study is to analyze the spatio-temporal changes in land cover and land use (1984-2010), as well to simulate future land cover for 2030 in the Neka River Basin (NRB), including the Hyrcanian forest, in northern Iran. For this purpose, we used detailed land cover maps for the years 1984, 2001 and 2010. The results showed that the highest deforestation occurred in the boundaries between forest and agriculture areas between 1984 and 2010. Comparing the observed and predicted land cover in 2010 yielded agreement of 96.41%. From 1984 to 2010, landscape metrics showed that the forest area evolved to more fragmentation, with less shape complexity and less connectivity. Projections for the future are consistent with observed changes for the Neka landscape, with a tendency to continue disaggregating and increasing diversity in a number of different patch types. Between 2010 and 2030, we observed the arrival of new crops, rangelands, and urban areas within the remaining areas of homogeneous forest. Changes in the Hyrcanian forest will cause alteration in ecosystem services, such as erosion control, water yield, timber harvest, and ground water reservation. Results of this work may represent a useful tool to provide strategies and territorial planning for sustainable management of the fragile Hyrcanian forest ecosystems in the Neka Basin. Modelización del cambio de la cobertura del suelo en los bosques de Hircania del norte de Irán: una perspectiva desde el análisis de la transformación y los patrones del paisaje RESUMEN. El objetivo principal de este trabajo es analizar los patrones espaciotemporales de cambios del uso del suelo (1984-2010) y generar escenarios para el horizonte temporal 2030 en la cuenca del río Neka, en el norte de Irán. Dicha DOI: http://doi.org/10.18172/cig.3279 © Universidad de La Rioja


Introduction
One consequence of world population growth is an alarming increase in land cover changes in the global environment, the conversion of rich ecological areas into cropland, and fertile lands into urban areas (Al-sharif and Pradhan, 2014). Land transformation, habitat degradation, and fragmentation are typical processes of landscape change (Paudel and Yuan, 2012). Pattern, rate, and trends of land cover change are vital to understanding forest dynamics, to enable conservation, and to evaluate management procedures (Gómez et al., 2011). Understanding trend and the amount of land cover change is an important issue among land use planners and environmental scientists, because these changes are one of the reasons for global environmental change, with deep ecological, hydrological, soil evolution, and erosion and social impacts. These changes are usually have an obvious anthropogenic source, but several topographic and geographic variables, such as slope, aspect, and altitude influence the nature and magnitude of these changes as well (Subedi et al., 2013).
Landscape matrices are important in understanding and monitoring anthropogenic effects on ecosystem dynamics. Landscape indices imply quantification of configuration, composition, shape, contagion, density, and diversity in patches of landscape. These indices divide into two main groups. The first is composition metrics, such as the number, density, and size of patches that correlate with the variety and abundance of patch types and nonspatial characterization of the landscape. Configuration indices attribute the spatial structure of landscape, such as neighborhood relationships and contagion (McGarigal and Marks, 1995;Rafiee et al., 2009). In the case of landscape ecology, configuration and composition landscape have important roles in determining natural habitats, hydrological processes, energy, and the nutrient cycle (Lee et al., 2009). Many studies showed an increase in the use of landscape metrics for monitoring of patterns of land cover changes (Lausch and Herzog, 2002;Antwi et al., 2008;Ramachandra et al., 2012;Plexida et al., 2014;Joorabian Shooshtari and Gholamalifard, 2015;Sakieh and Salmanmahiny, 2016;García et al., 2017). Boentje and Blinnikov (2007) assessed forest loss in the Green Belt around Moscow, Russia during the 1991 to 2001 period. They found that conversion of forest to suburban residential and commercial comprised about 14.6% of the land within 20 km from the beltline. The results of landscape indices showed that forest patches became: 1) larger on average during the 1991-2001 time interval; 2) a reduction in forest edge observed in 2001 compared to 1991; and 3) less connectivity between forest patches. Xin et al. (2014) studied land use and landscape change pattern using remote sensing and landscape metrics in Turkmenistan between 1976 and 2011. The metrics proposed increases in landscape fragmentation and declining connectivity of patches. Kamyab and Salman Mahini (2014) analyzed spatio-temporal patterns of urban changes in the city of Gorgan, Iran. Between 1987 and 2001, the city of Gorgan displayed an increase in size, simpler in patches shape, and less fractal dimension.
In addition to landscape metrics, indices of spatial processes for landscape transformation provide insight and comprehensive information about changes in landscape. Perforation, dissection, fragmentation, shrinkage, and attrition relate to a decline in area and involve land cover degradation. Aggregation, creation, and enlargement are related to an increase in area and creation of new land cover units. Shift and deformation do not relate to area change (Bogaert et al., 2004;Vranken et al., 2011). Munsi et al. (2010) evaluated land cover change in the Dehradun forest, Uttarakhand, India, using landscape metrics and landscape transformation type. The transformation type was dissection of dense forest, and creation of open and scrub forest during the periods 1990-2001 and 2001-2006. The spatial process of landscape transformation in opencast mining lands was explored by Fagiewicz (2014). Analysis of forest, settlement area, and arable lands indicated that the landscape transformation process was fragmentation, that natural reservoirs represented shrinkage, and that wetlands showed attrition. Redondo-Vega et al. (2017) explored land cover changes caused by mining activity in the northwestern mountains of Spain, using 1956-57 aerial photographs and 2014 colour orthophotos. Mining operations caused a radical change in cultivated surfaces, livestock grazing, and low hillsides. Additionally, the mining of mineral resources altered or eliminated valuable features of geomorphologic heritage.
In addition to landscape indices and change process, land cover change models can be useful tools for understanding dynamic future landscape. GIS models and RS data can be used to simulate how landscape elements changes over time and place and to explore various types of scenarios for the future (Paudel and Yuan, 2012). Hyrcanian forests in northern Iran are a source of biodiversity, diverse habitats, and commercial timber production, and have a vital role in soil conservation, carbon storage, gentler air, and water purification (Poorzady and Bakhtiari, 2009;Joorabian Shooshtari et al., 2012). Transition in these forests to another land cover type is one of the major challenges of recent years; therefore, their protection is important and necessary (Kelarestaghi and Jafarian Jeloudar, 2011). Losses in these forests have major effects on biodiversity and ecosystem services. Consequently, the analysis and prediction of Hyrcanian forest cover dynamics during different periods is important to improve comprehensive understanding about the quantity of deforestation in these valuable ecosystems. The main goal of this study is analyzing landscape transformation type from 1984 to 2030 and performance assessment using SimWeight (similarity weighted instance-based learning algorithm) to discern modeling change potential in the Neka River Basin (NRB), northern Iran.

Study area
The study area is in Iran's Neka River Basin, which consists of 1871 km 2 in Mazandaran and Golestan Province (Fig. 1). The elevation of the area ranges from 50 to 3791 m a.s.l., with average elevation of 1576 m a.s.l. The mean annual precipitation is 600 mm and the mean annual temperature is 17ºC (Ghanbarpour et al., 2014). Almost 44% of the study area is covered by rangelands. Hyrcanian forests occupy about 36% of the Basin, followed by agriculture area occupying over 18% in 2010. The Hyrcanian forests zone covers the Talysh Mountains in southeastern Azerbaijan and the northern slope of the Alborz Mountains in Iran. These forests are composed mostly of deciduous forest with high biological diversity and endemic species. In total, 3234 species belonging to 856 genera and 148 families of vascular plants have been reported in these forests (Mostajeran et al., 2016). Human activities have severely reduced the areas of Hyrcanian forest from 3,600,000 ha in 1942 (Saei, 1942) to 1,920,000 ha in 1990 (Moshtagh Kahnamuii and Rasaneh, 1990). The study area is potentially at high risk of flooding (Mohammadi et al., 2014), as occurred in 1999.

Historical changes, transition potential, and future modeling
The input data required for Land Change Modeler (LCM) includes at least two land cover maps at different time periods (Eastman, 2012). In the present study, in order to explore changes analysis, we have used the land cover maps from 1984, 2001 and 2010, which had been extracted from satellite images. The amount of conversion, spatial distribution of transitions, gains and losses bet ween land cover categories were calculated and analyzed between the periods 1984-2001 and 2001-2010 by LCM for Ecological Sustainability (Schulz et al., 2010). We conducted land cover change model development through the following steps. First, we used the land cover maps of 1984 and 2001 in developing a land cover map for the year 2010. The actual 2010 land cover map was applied to evaluate the performance of the model based on the Kappa Index of Agreement (KIA) (Pontius et al., 2001). This stage was adopted in order to evaluate the model simulation. Then, land cover maps of the years 2001 and 2010 were performed to create prospective land cover maps of 2030. SimWeight was used to develop suitability maps of the likelihood for transition from one class to another (Sangermano et al., 2012;Lin et al., 2014). It uses a slightly modified variant of the K-nearest neighbor procedure. Sample size and K parameters must be specified (Eastman, 2012). In SimWeight, each variable's importance in each submodel is computed using the relevance weight indicator (Eq. 1), where SD is standard deviation of the variable (Sangermano et al., 2010). We ignored variables with an estimated relevance less than 0.01 in modeling each submodel.
Performance of SimWeight was evaluated using Peirce Skill Score, which is the difference between a hit rate (calculates the mean transition potential among the validation pixels that actually went through the transition) and a false alarm rate (calculates the mean transition potential among pixels that were eligible to change but did not) (Eastman, 2012). An extensive description of SimWeight can be found in Sangermano et al. (2010). LCM identifies the transition potential maps based on historical changes and explanatory (static or dynamic) variables (Lin et al., 2014). We selected 9 variables, including elevation, slope, empirical likelihood to change map (which is the empirical probability of change between 1984 and 2001), distance to residential areas, forest areas, agricultural lands, rangeland, major road, and fluvial streams (Fig. 2). Distance maps were prepared with values indicating Eucidian distances from target layers. The 30 m digital evaluation model (DEM) from ASTER Global DEM (METI, Japan and NASA, USA) was used to prepare the slope. Cramer's V coefficient (based on the chi-square statistic) was examined in order to explore the explanatory power of these variables (Pérez-Vega et al., 2012). A Cramer's V value of about 0.15 or higher, and 0.4 or higher, reveals that the potential of variables is useful and good, respectively (Eastman, 2012). Six types of transition were produced, including 1) transition from forest to agriculture area; 2) forest to residential area; 3) forest to rangeland; 4) agriculture to residential area; 5) agriculture to rangeland; and 6) rangeland to residential areas. The quantities of change from one class to another class were calculated by Markov chain analysis (Khoi and Murayama, 2010). After calculating future demand for land, multi objective land allocation (hard prediction) was run to generate the future land cover (in this case, for the years 2010 and 2030).

Calculation of landscape metrics and landscape transformation analysis
Seven indices were selected to analyze the landscape transformations: Patch Density (PD), Perimeter Area Weighted Mean Ratio (PARA_AM), Interspersion and Juxtaposition Index (IJI), Patch Cohesion Index (COHESION), and Mean Euclidian Nearest-Neighbor Distance (ENN_MN) were analyzed at the class level, and two metrics including Shannon's diversity index and Contagion were calculated at the landscape level using the FRAGSTATS program (UMass Landscape Ecology Lab, Amherst, USA) ( Table 1) (McGarigal et al., 2002). PARA rises with increasing patch shape complexity. The IJI index represents the degree to which patch types are interspersed. Contagion index measures the degree of land cover types clumping in a landscape (Riitters et al., 1996). The value of this index is 0, when the patch types are maximally disaggregated, and 100 when all patch types are maximally aggregated (McGarigal et al., 2002). A decision tree algorithm within LCM was applied to identify the spatial process in landscape transformation (Fig. 3) (Bogaert et al., 2004). It consists of number of patches, perimeter, and area of the class in two different time periods.

Analysis of land cover change
Analysis of the land cover changes in the NRB from 1984 to 2010 are depicted in Table 2. The area of different classes indicates that the dominant land cover types are rangeland and forest, representing 80% of the overall case study.  years 1984, 2001 and 2010 and area changed during the 1984-2001, 2001-2010 and 1984-2010  The analysis revealed that major decline, 3000 ha, occurred in forest (1.6% of the total area and 4.1% of the forest area) and 1971 ha (1.1% of the total area and 2.8% of the forest area), respectively, between the periods 1984-2001 and 2001-2010. During the period 1984-2001, forest converted to agricultural lands and rangelands by 2260 ha and 601 ha, respectively. The GIS analysis reported that built-up area and road expansion in the NRB was 94 ha (58.8% of the residential area) and 189 ha (9.3% of the road area), respectively, in the period from 1984 to 2001. Transition from agricultural land to residential area mainly contributed to the increase in this class. During the second time period, 2001 to 2010, the major changes in the study area occurred in forest, with conversion to agricultural areas (1018 ha) and to rangelands (885 ha). Another change was the increase in built-up area and road between 2001 and 2010. The study revealed that the area occupied by water bodies and bare lands almost remained constant during the 26-year period from 1984 to 2010. Spatial patterns of changes from forest to agricultural areas and rangelands from 1984 to 2010 are shown in Figure 4.

Prediction of land cover changes
The results of Cramer's V statistic, representing the association level of the explanatory variables and the existent land cover in the NRB, are shown in Table S1  (see Supplementary Material, Table S1). Empirical likelihood to change and slope are the variables that showed the highest and lowest relationship with land cover changes. Transition potential modeling was carried out using SimWeight for 6 submodels in the study zone. Results of the relevance weight of independent variables and Peirce Skill Score for SimWeight validation are presented in Tables 3 and 4, respectively. Nine variables showed different levels of importance in each submodel. All submodels were strongly affected by the empirical likelihood to change variable. The lowest Peirce Skill Score (0.37) was for conversion from agriculture areas to residential submodel and the largest value (0.84) was for transition from the rangelands to residential areas submodel. The transition probability matrix calculated by Markov chain analysis in the NRB is presented in Table S2 (see Supplementary Material, Table S2).  The highest likelihood of transition from forest to agriculture occurred between 2010 and 2030. The contribution to increasing residential areas comes from agricultural lands. The 2030 land cover map showed a reduction in forest (4225 ha) and increases in rangeland (2024 ha), agricultural areas (2144 ha), and residential (57 ha), as compared with 2010 land cover (Fig. 5). The predictive performance of the LCM was validated by comparing the model simulating 2010 land cover with the observed 2010 land cover; agreement was 96.41% (Fig. 5).

Analysis of landscape metrics and spatial processes of landscape transformation
Patch density (PD) in forest class increased from 1.5 in 1984 to 1.9 in 2030. This index for forest recorded a decline during all time intervals except 2010-2030. The number of patches per unit area and the reduction in area in the forest class between 1984 and 2010 declined, attributed to attrition in this class. During the 2010-2030 period, forest represents dissection because of an increase of patch density (from 1.3 to 1.9) and a decrease in area (from 67198 ha to 62973 ha). Changes in PD (from 1.62 to 1.51) showed decreasing and increasing areas under agricultural lands, referring to the aggregation transformation type for the 1984-2010 period (Fig.  6). Rangelands showed results similar to those of agricultural areas, with a spatial process of aggregation and creation during 1984-2010 and 2010-2030, respectively. The number of patches per unit area in residential use declined between 1984-2001, while PD increased between 2001and 2030 in addition, the increasing area in this class means aggregation for the first period and a creation transformation process for the second and third periods (Table 5). PD values of road and area under this class confirmed an increase, which reflects that the change process was creation from 1984 to 2010. The change process type for road was dissection between 1984 and 2010. Residential areas underwent creation, indicating that PD (from 0.10 to 0.15), and area (from 253 ha to 334 ha) increased (Table 5). PARA is a simple measure of shape complexity, but without standardization to a simple Euclidean shape (McGarigal et al., 2002). Rangeland and road for PARA showed a decreasing trend during the 1984-2030 period, which represents less shape complexity. Residential, agriculture areas, and forest decreased during the 1984-2010 period and increased between 2010 and 2030, indicating more shape complexity of these categories in the third period (Fig. 6). The COHESION metric reflects the physical connectedness of the particular land use class (Lee et al., 2009). Connectivity is important to the constancy of both biotic populations in a fragmented ecosystem (Joorabian Shooshtari and Gholamalifard, 2015). Between 1984 and 2030, this index increased for agricultural lands, residential, rangeland and road, thus making them more connected. Connectivity between patches of forest from 2010 to 2030 revealed reduction (Fig. 6). Analyses of COHESION metric deduced a decline in forest from 1984 to 2030, which showed less physical connectedness between patches in the NRB. High values of IJI indicate patch types that are equally adjacent to each other (McGarigal et al., 2002;Munsi et al., 2010;Lin et al., 2007). Decreasing of this metric for forest and residential areas becomes more clumped over the years 1984-2030, while for agriculture areas and rangelands, it increases (Fig. 6). The lowest value of the mean Euclidian nearest-neighbor in the NRB was obtained in road class. This class displayed a decreasing trend (42.2 to 41.8) from 1984 to 2030. Values of ENN in forest decreased between 2010 and 2030, indicating that forest patches became more isolated, affecting the movement and dispersal of species (Midha and Mathur, 2010). A decreasing trend in ENN for the agriculture area (88.8 to 80.8) and residential (109.4 to 80.9) was apparent from 1984 to 2030. A decrease in the score of ENN was observed for rangeland from 80.1 in 1984 to 85.1 in 2030. Moreover, by analyzing landscape metrics at the class level, assessment of patterns of the NRB was conducted at the landscape level. The CONTAG metric for the entire NRB landscape showed reduction from 66.9% in 1984 to 66.5% in 2030. The increase in SHDI index was seen for the landscape between 1984 and 2030 (Table 6). These changes showed the NRB landscape's tendency to disaggregate with increasing diversity. According to very small changes in bare land and water body classes, the focus is not on their results, but upon graphic display (Fig. 6).

Discussion
In the present study, GIS modeling and landscape transformation analysis were innovative planning tools to explore land cover dynamics and spatial change process of NRB landscape transformation from 1984 to 2030. The results of change analysis over the 1984-2010 study period demonstrated forest loss of 4970 ha. Logically, the highest deforestation occurred in the surrounding agricultural lands, because of their accessibility. Additionally, the increase in population growth, use of wood as fuel, construction, and repairing villagers' houses are other reasons for deforestation in the study area (Joorabian Shooshtari et al., 2010). Abdullah and Nakagoshi (2006) reported agricultural development as a major variable in landscape dynamics. The spatial expansion in rangelands (1502 ha) and agriculture areas (3053 ha) was observed between 1984 and 2010. A slight increase in the extent of residential areas and transportation was revealed during the 26-year period from 1984 to 2010. The increase of rangeland was mainly due to the reduction of forest during the 1984-2010 period.
We applied SimWeight in order to conduct transition potential modeling in 6 submodels. In our study, the highest association was observed between empirical likelihood to change map and elevation, and land cover changes based on Cramer's V coefficient. Munsi et al. (2012) reported that the empirical likelihood of change and elevation showed strong correlation with land cover changes for modeling patterns of forest cover in India. Bagheri and Shataee (2010) performed modeling of forest dynamics in the Chehl-Chay basin, Golestan Province, Iran. Their results showed that elevation increase causes an increase in deforestation, because of more villages and road development in higher altitudes. In the NRB, proximity to forest, residential areas, rangeland, and agriculture variables showed overall Cramer's V > 0.23, meaning that these independent variables are especially important drivers in the description of land cover change dynamics. Distance to residential and road variables revealed high relevance weight in the forest to residential submodel and distance to rangeland was highly important to the transition from forest to rangeland. Eraso et al. (2013) verified that proximity to road, cities, and pastures were the main factors affecting deforestation of montane forest. Slope showed high values of relevance weight in the transition from forest, agriculture, and rangeland to residential submodels, as steep areas are not in expected to support more anthropogenic activities, and are difficult to cultivate or urbanize in this area; therefore, this variable is more important in these submodels. Distance to river affects the location of agricultural areas; accordingly, the proximity to water bodies is related to deforestation (Khoi and Murayama, 2010).
In the NRB, the results of SimWeight validation showed reasonable performance. Sangermano et al. (2010) mentioned that SimWeight is an effective procedure for calculating transition potential, without requiring specification of complex parameters, and a distribution-free non-parametric method. In the present study, a Markov chain matrix was applied to determine the quantity of changes from each to another category for future dates (Kim, 2010;Mas et al., 2014;Camacho Olmedo et al., 2015). Forest showed the highest probability of changes to agriculture and rangelands based on Markov transition probabilities from 2010 to 2030. The KIA indicated a good agreement between the 2010 model's predicted land cover and 2010's actual land cover as a reference in the NRB. According to our prediction results, deforestation will probably take place in buffer areas surrounding agriculture areas and rangelands by the year 2030. A slight increase was observed in residential areas between 2010 and 2030.
Moreover, in the present study, analysis of landscape pattern was conducted by landscape metrics and landscape transformation analysis. While numbers of patches have increased, the area under the class has been reduced, indicating that the transformation type is dissection or fragmentation. If the number of patches has decreased while the area of the class remains constant or increases, this indicates that the transformation type is aggregation (Munsi et al., 2010;Raja Naqvi et al., 2014). In the present study, during the 1984-2001 and 2001-2010 periods, transformation type was attrition for forest because of increasing human pressure. Between 2010 and 2030, forest was dissected, which means that the number of patches increased while the area under the class (67198 ha to 62973 ha) is decreasing. Residential areas showed aggregation during 1984-2001, representing more clumping and compacting in this class. Population growth causes increasing demand for agriculture fields and built-up area, which leads in the NRB to the decline of forest and ecosystem degradation. Landscape transformation type in agriculture areas was aggregation and creation between 1984-2010 and 2010-2030, mainly because of removal of forest patches. Additionally, investigation of results obtained from landscape metrics showed that increasing the SHDI index combined with a reduction of CONTAG revealed that scenery becomes disaggregated and increasingly diverse from 1984 to 2010. This tendency is expected to continue over the next decades.
Hyrcanian forests are unique, very old with a very long history and some relic species date back to the Tertiary period (Rouhi-Moghaddam et al., 2007), so slight and limited changes in these forests are very important and have consequences for ecological condition, basin hydrology, and sediment transport. The NRB is prone to flooding, which could contribute to the decline in the forest area, followed by a decrease in rainfall retention in soil and increased runoff (Talebi Amiri et al., 2009). According to the slight change in land cover of the Neka Basin, there is apparently no major driver of hydrological changes in the region, except perhaps for the increased water demand for agricultural areas and the water supply for residential and industrial regions. Sometimes, a very small change in land cover and landscape structure can have a significant impact on regional soil erosion and sediment yield (García-Ruiz et al., 2008). The recommendation for future studies is explore the response of sediment transport to land cover change in the Neka Basin. Furthermore, many channels extracted the water of Neka River for agricultural use; it is important to obtain comprehensive information about the increase sediment yield within channels, in order to make informed management decisions.

Conclusions
The results of this study suggested that landscape indices and landscape transformation analysis are effective tools for understanding the ecological patterns of land cover change and maintaining and managing the natural resources in the study zone. In general, forest in the NRB showed dissection transformation type during the 2010-2030 period. Deforestation in the present case study has effects on the reduction in essential nutrients, regulation of runoff and regional climate, rainfall infiltration, and soil conservation (Joorabian Shooshtari and Gholamalifard, 2015). It is necessary to develop strategies for landscape composition, configuration, structure, and restoration for the Hyrcanian forest in the NRB. Table S1. Results of Cramer's V statistic, revealing the degree of association between the explanatory variables with distribution of land cover classes for LCM.