Simulating Spatiotemporal Changes in Land Use and Land Cover of the North-Western Himalayan Region Using Markov Chain Analysis

: Spatial variabilities and drivers of land use and land cover (LULC) change over time and are crucial for determining the region’s economic viability and ecological functionality. The North-Western Himalayan (NWH) regions have witnessed drastic changes in LULC over the last 50 years, as a result of which their ecological diversity has been under signiﬁcant threat. There is a need to understand how LULC change has taken place so that appropriate conservation measures can be taken well in advance to understand the implications of the current trends of changing LULC. This study has been carried out in the Baramulla district of the North-Western Himalayas to assess its current and future LULC changes and determine the drivers responsible for future policy decisions. Using Landsat 2000, 2010, and 2020 satellite imagery, we performed LULC classiﬁcation of the study area using the maximum likelihood supervised classiﬁcation. The land-use transition matrix, Markov chain model, and CA-Markov model were used to determine the spatial patterns and temporal variation of LULC for 2030. The CA-Markov model was ﬁrst used to predict the land cover for 2020, which was then veriﬁed by the actual land cover of 2020 (Kappa coefﬁcient of 0.81) for the model’s validation. After calibration and validation of the model, LULC was predicted for the year 2030. Between the years 2000 and 2020, it was found that horticulture, urbanization, and built-up areas increased, while snow cover, forest cover, agricultural land, and water bodies all decreased. The signiﬁcant drivers of LULC changes were economic compulsions, climate variability, and increased human population. The analysis ﬁnding of the study highlighted that technical, ﬁnancial, policy, or legislative initiatives are required to restore fragile NWH regions experiencing comparable consequences.

Land 2022, 11, 2276 3 of 18 our knowledge, no prior research has applied the CA-Markow model to forecasting LULC in the region under study. The aims of the current research are threefold: (1) to examine the spatial and temporal shift in LULC in the Baramulla area from 2000 to 2020; (2) to use the CA-Markow model to project LULC for the year 2030, and (3) to demonstrate how this movement will impact future land-use planning. To predict changes in landscape structure and use natural resources responsibly, regulators, resource managers, and government bodies should keep a careful eye on these spatio-temporal situation models.

Study Area
Baramulla is located in the Northern-Western part of the Himalayan Mountains between 34 • 16 and 34 • 40 N latitude and between 73 • 45 and 75 • 35 E longitudes. The study area's climate is temperate, with moderate-heavy snowfall. The average temperature of the study area is 25 • C, and the annual average rainfall ranges from 860 to 1050 mm. The highest precipitation is received from March to July, with March being the wettest month. The study area (7000 Km 2 ) is in the northmost part of India, sharing a boundary with Pakistan, and is located at an altitude of 1060 to 5200 m.a.s.l. The altitude on the north side is lower than in the center and other parts [45]. The topography is complex and constitutes both valleys and mountains, with most of the area under the mountain landform. The primary source of irrigation water is the Jhelum River and its tributary, and on its northeast side is Wular Lake, which is Asia's biggest freshwater lake. The primary soil type of the study area is silty clay loam and is grouped into alluvial and forest soil. The main soil orders are Alfisols, Inceptisols, and Entisols [46,47]. The vegetation types in the study area include natural forests such as Pinus spp., Ciderus spp., Salix, Populus, and Fir spp. (Figure 1, Table 1).

Datasets
In the current study, multi-temporal satellite data were used in this research to model and predict LULC. Landsat 7 Enhanced ETM+ data dated 5 May 2000 and 26 May 2010; and Landsat 8 OLI data dated 5 May 2020 were downloaded from the website of USGS (www.earthexplorer.usgs.gov, accessed on 15 October 2022). Landsat images having a resolution of 30 m were used for the years 2000, 2010, and 2020. The significant steps used in the image analysis included (a) data collection, (b) classification of LC during these years, and (c) prediction of LC using various models. The advanced space-borne thermal emission and reflection radiometer (ASTER) DEM having a resolution of 30 m, was downloaded from (https://terra.nasa.gov/ accessed on 31 January 2022) for ortho-rectification [48] [ Altaf et al., 2013]. The maximum likelihood classification approach was used to differentiate various classes, which included agriculture land (paddy, maize), horticulture (apple, pear, cherry, walnut), forest LC (pinus, cederus, populous, and salix), construction (houses, roads, schools), fallow lands, water bodies, snow, urban, and vegetation. ERDAS imagine (v. 9) was used for supervised classification. We used ArcGIS 10.2 to assess the changes in LC using the following formula:   A (%) = percentage of difference in land use and land cover type between the initial and final value, At 1 is the initial area of land use and land cover, and At 2 is the area of land cover and land cover at the final stage. The rate of land use and land cover changes was calculated by the following formula: where; ∆R = rate of change, C = current area coverage of LC, P = previous area coverage of LC, T = time interval between C and P Land 2022, 11, 2276 5 of 18

Supervised Classification
The Landsat images were classified into different classes based on the requirements of the study area. Eight major LULC classes, viz., agriculture, horticulture, forest, urban, waterbody, fallow land, other vegetation, and snow were identified for mapping the area ( Table 2). Maximum likelihood classification was used to map all the land use/cover classes [49,50]. Prior to choosing the training samples, a thorough empirical investigation of the local toposheet, Google Earth imagery, and satellite data was conducted. The minimum number of training samples for most classes was 100. With the use of GPS and local guides, a field survey was carried out in several locations, encompassing all of the LULC classes, to verify the correctness of any questionable places on the ground. Due to mountainous topography, rough terrain, and steep slopes, few areas were not accessible [51].

Methods
Changes in LULC are related to the socioeconomic and biophysical factors existing at a specific time period and place. The LULC changes in a particular moment drive the LULUC change at a certain potential area, and thus modeling takes that scale of change into account. Any shift in LULC is based on a combination of socioeconomic factors and biophysical features. Biophysical parameters include temperature, rainfall, radiation, water availability, soil depth, soil fertility, and nutritional status, whereas socioeconomic conditions include the farmer's financial situation, people's preferences, and tastes [48]. LULC changes at the global level are due to variations of some macro-variables like technology level, quality of life, and population, which influence environmental change in the long run. Still, at the local level, these changes are due to factors such as changing economic structure of farmers, topographical conditions, relative desirability aspect, farmers' preferences for a specific crop, migratory tendencies, changing climate, and bioresource availability [52]. The following are some of the major LULC change models that are worth discussing in this paper.

Cellular Automata
Ulman, in 1940, developed the most dynamic LULC model known as cellular automata (CA), and Tobler was the first to determine its use in LULC and the modeling of geographic data [53,54]. A CA model is a spatio-temporal model for analyzing future LC. This model is distinct in time and space and carries out complex space-time simulations. The CA model consists of raster space, grid cells, neighbors, time, and transition rules [7,55]. The grid cell and its neighbors determine the set of states and neighborhood arrangement of the cells. The transition rule analyzes the state transition of each cell as a function of the position of its neighboring cell. The time component then updates the configuration and composition of each cell simultaneously. In a CA model, the state L c+1 for every cell is decided by the cell itself and by the state Lc of neighboring cells, implying that the transition rule analyzes changes. where L is the set of states of each cell, and N is the neighboring cells. The c and c + 1 are different moments, and f is the rule of the transformation.

Markov Prediction Model
Andrei A. Markov, a Russian mathematician, gave the Markov chain analysis model in 1970, but Burnham's use in agricultural science was first imposed in 1973. Markov's model is speculative, uses patterns that provide change among various land exploits, and depicts future trends. A Markov model is a complex operation that is distinct in time and space. This model provides LULC images, a transition probability matrix, and a transition area matrix and predicts LULC change trends [56,57]. The transition probability matrix best explains the LULC change over time and space. The transition area matrix determines the number of pixels predicted to change in each LULC type over a specified time unit. The most important application of the Markov model is to determine the transitional pattern between various LULCs by observing the transition probability between the primary and terminal states. The Markov model is a raster data model of GIS, which predicts uninterrupted data over space. In this model, the area is split into pixels, and each pixel is filled with measured values. The cell values are written in the columns. The pixel values represent land use and land cover change data from satellite images [58,59].
The Markov model is a raster data model of GIS, which predicts uninterrupted data over space. In this model, the area is split into pixels and grid cells, and each pixel or grid cell is filled with measured values. The cell values were written in the columns. The pixel values represent LULC and LULC change data derived from satellite images. The Markov model could be represented by the set of states L = (L0, L1, L2, L3 . . . .. Ln), predicting that the current state is Lc which gets changed to state Ld at the adjacent step with the transition probability denoted by Pij [60,61]. The state Lc+1 in the system could be obtained by the previous step Lc using the formula [62].
where n is the LC type number, L is LC type.

Land Cover Change Detection
The change in the area is determined by the change in the image grid cell with a unit size of 30 × 30 m. The cell space represents the whole land-use area. The IDRISI 17.0 Markov model was used to obtain the state transition probability matrix and LULC transfer matrix. The multi-criteria evaluation (MCE) suitability maps and multi-objective decisionmaking module of IRDISI were used to obtain suitability. The results were derived under a 5 × 5 pixels filter, meaning that the central cell is affected by its 5 × 5 m neighboring cells. The LC prediction for 2030 was then carried out using the Markov modeling IRDISI 17.0.

Model Validation
Detection and validation are the core steps in the modeling process, especially in detecting future changes where no datasheets are available. The output of a model validation determines the model's efficiency ( Figure 2). The generated 2000 and 2010 land use cover was used to simulate the 2020 LULC using the CA Markov model. The simulated 2020 LULC was validated using the actual 2020 LULC. This output value is acquired by comparing the results of the simulation with the present-day changes using kappa statistics, e.g., kappa for quantity (Kquantity), kappa for location (Klocation), and kappa for no infor-  [58,63]. After acceptable results are obtained from the model using kappa statistics, the model is considered fit for predicting future changes. The Kno represents the standard kappa index's variation and determines the simulation run's overall accuracy. The Kquantity and Klocation indicate the quantity and location of the simulation. The value 1 of these three indices represents the perfect simulation, and the value 0 represents an imperfect simulation. The VALIDATE module was used to validate the predicted data with reference data. Cross-tabulation was used to determine the predicted data with actual data to obtain the Kappa index.

Dynamic Rate Changer of Land Use and Land Cover
The land cover dynamic degree represents the severity of land-use change in a given period, which forms the important index of land-use change (Mukherjee et al. 2009). This is measured by the single and comprehensive LC dynamic degree. A dynamic degree spatial analysis model of LULC was also used [62]. is measured by the single and comprehensive LC dynamic degree. A dynamic degree spatial analysis model of LULC was also used [62].
where L i and L are the dynamic degrees of single and comprehensive LULC dynamic degrees, LA (i, c 1 ) is the land-use area of a particular LC type at an initial stage, and LA (I, c 2 ) is the land-use area of a particular LULC type at a later stage, ULA i is the area not changed, c 1 and c 2 represent the initial and final time. ToR is a transfer-out rate, and TiR is a transfer-in rate, LoR is a sum of ToR and TiR.

Image Classification and Accuracy Assessment
Classification results of the Landsat images in 2000 and 2020 were used to determine the potential area of LC change. Figure 3 represents the classification results of Landsat images in 2000, 2010, and 2020. Table 2 shows the LULC statistics for the three periods. These tables and figures show that the study area has been classified into eight major LULC classes: agriculture, horticulture, forest, fallow, snow, vegetation, water bodies, and construction. After the spectral classification, the area of each class and their area percentage were measured. The results showed that forests covered the most significant area across the years, with an area proportion of 26.15%, 24.50%, and 24.31% in 2000, 2010, and 2020, respectively. Agriculture was the second most dominant LC, with an area coverage of 22.30, 21.29, and 20.04%, with a declining pattern in 2000, 2010, and 2020, respectively. The study area is known as the horticultural belt of the North-Western Himalayas, and the proportion of area coverage increased gradually from 3.92 to 5.77% between 2000 and 2010.
The proportion of construction area increased at an incredible pace from 2000 to 2020, with a percentage of 9.71% in 2000, 13.17% in 2010, and 14.11% in 2020, showing the urbanization and improvement in the economic condition of the inhabitants. The water bodies are degraded at the expense of construction and other farm activities. The area of water bodies in 2000 was 1.06%; in 2010, it was 0.98%; and in 2020, it was 0.77%. The width of rivers decreased sharply. The data for validation was manually and randomly selected based on Google maps. Table 2 pertains to the results of processed images of three periods. Both the user's accuracy and the producer's accuracy are obtained from the confusion matrix. The overall percentage of accuracy was 83.12 in 2000, 90.37 in 2010, and 90.87 in 2020 (Table 3).
across the years, with an area proportion of 26.15%, 24.50%, and 24.31% in 2000, 2010, and 2020, respectively. Agriculture was the second most dominant LC, with an area coverage of 22.30, 21.29, and 20.04%, with a declining pattern in 2000, 2010, and 2020, respectively. The study area is known as the horticultural belt of the North-Western Himalayas, and the proportion of area coverage increased gradually from 3.92 to 5.77% between 2000 and 2010.     Tables 5 and 6 reflect the dynamic change in these three periods. From 2000 to 2010, the greatest net change was observed in horticulture, followed by urban and fallow land, with change percentages of 47.25%, 35.59%, and 19.10%, respectively. The amount of the urban and fallow land-use area transformed into another LULC during 2000-2010 was 241 and 170 km 2 . The amount of horticultural land-use area converted from another LC was 129 km 2 . The lowest percent area converted into another was vegetation, followed by water bodies and snow, with area proportions of 2.73, 8.11, and 10.75, respectively. Due to the harsh climate, the snow area is difficult to change into other land-use types. The amount of forestland transformed into other LULC was 115.00 km 2 , accounting for 6.31%. From 2010 to 2020, the horticultural land cover had the greatest change, with an area of 186 km 2 and a proportion of 46.27%. The second highest change was observed in the water class, followed by fallow land with area coverage of 14.00 and 143.00 km 2 , having area proportions of 20.59 and 19.86%, respectively. As a result of minimal snowfall and ongoing sand excavation, the area's water covering shifted throughout this period. From 2000 to 2010, the percentage change in horticultural land cover was the biggest, at 53.57 percent, with an area of 3.15 km 2 . Other land cover types that changed included water, vegetation, snow, forest, urban, and agricultural. The next two largest land uses to experience a shift were fallow land (54.2%) and urban land (31.2%), each accounting for 313 km 2 and 307 km 2 , respectively. Horticultural land was transformed from fallow land, followed by urban, water, snow, forest, agriculture, and vegetation. Tables 5 and 6 show the dynamic shifts in LULC from 2000 to 2010 and 2010 to 2020, respectively. With transfer areas of 235.28 and 205.00 kma −1 and transfer rates of 2.64 and 2.85 percent, respectively, fallow land experienced the biggest rate shift. The second highest transformation was in the water area with 2.14 and 2.83 km a −1 from 2000 to 2010 and 2010 to 2020. This was followed by snow LC with a transfer rate of 1.74 and 0.97 kma −1 . The lowest transformation rate was observed in vegetation with a value of 0.35 and 0.12 kma −1 . The pattern of transformation rate was Fallow land > Water > Snow> Urban > Horticultural > Forest > Agriculture > Vegetation. The gain rate was highest in horticultural LULC from 2000 to 2010 and 2010 to 2020, with 4.65 and 2.72 kma −1 . The lowest gain rate occurred in a forest area with a value of 0.6 and 0.19 kma −1 . The pattern of gain rate was Horticultural > Urban > Water > Fallow land > Snow > Vegetation > Agricultural > Forest. The dynamic degree followed the pattern of Water > Snow> Fallow land > Horticultural > Agriculture > Forest > Urban > Vegetation. The dynamic degree is meager than the rate of change. The data depicts intensive LC change due to gain, transfer, dynamic, and change rates. A low value of dynamic rate than the rate of change depicts that no single land is predominant. Tables 7 and 8 represents the land use status between the years 2000 and 2010. We can quantify the varying transfer probabilities using Markov's transformation matrix, elucidating several LULC transfer procedures. The rows represent the land-use change rate during 2000-2010, while the columns represent the land-use status and transfer during that period. The highest gain in the area was found in horticulture LC, whose net area was 273 km 2 , which accounts for 3.92%. The net increase in horticultural land is due to forest land transfer-in, agricultural land transfer-in, and fallow land transfer-in with transition probability of 0.1859, 0.2136, and 0.1479, respectively. The second-highest increase was in urban land, having an area of 677 km 2 , which accounts for 9.71. Its expansion is due to fallow land transfer-in and agriculture land transfer-in having a state transition probability of 0.1430 and 0.0159, respectively. The third highest increase was in vegetation LC, and its enhancement was due to water area transfer-in and snow area transfer-having a state transition probability of 0.0196 and 0.0476. Water and snow area decrease. The water area decreased by 16 km 2 , and the snow area decreased by 43 km 2 , accounting for 8.11 and 10.75%. The water and snow area are transformed into vegetation and forest LC. Tables 9 and 10 represent the LC status and probability matrix from the year 2010-2020. The highest gain in the area was found in horticulture, whose net area was 402 km 2 , which accounts for 5.77%, having a transition probability of 0.2391. The expansion or gain is due to agriculture and forest land transfer-in. The second-highest increase was in urban LC, having an area of 918.00 km 2 , which accounts for 13.77%, and its expansion is due to agriculture horticulture land transfer-in, agriculture land transfer-in, vegetation land  We generated the transition area and probability matrix in CA-Markov model IDRISI software. By running the CROSSTAB module in IDRISI, the actual LC map of 2020 is the foundation for a prediction map. Using the kappa index, we evaluated how reliable the findings were. As the Kappa index goes from 0 to 1, with a value of 0.5 indicating a lack of connection between the predicted and correctly classified photos, we can say that the predictive images have a low level of accuracy. When the kappa index is between 0.5 and 0.75, significant shifts and correlations can be seen between the two images. We observed an excellent correlation between the two sets of data (predicted and observed maps) since the Kappa index in our study was greater than 0.75. The estimated LC map for 2030 was based on a greater degree of kappa values, a transition probability matrix, and a transition area matrix (Table 11).

Historical Land Use-Land Cover Changes
A spatial data analysis approach was used to understand the spatio-temporal change in LC on decadal time series in the North-Western Himalayas. A Significant difference was found in North-Western Himalaya's multi-temporal LC map from 2000 to 2020. The results of image classification and the patterns of land cover changes are shown in Figure 3. There has been a considerable shift in the area between LULC types in North-Western Himalayas, as seen in Tables 3 and 5. Researchers have observed similar land cover change patterns across different regions [64][65][66]. This research found that changes in agricultural land were related to urbanization and horticultural land use demand.
Agriculturally productive land decreased by 20.02 percent between the years of 2000 and 2020. The expansion of horticultural and built-up classes was related to a decrease in snow, vegetation, and agricultural land use. These results have been reported by various researchers working in the Himalayas. These results align with several studies in which agricultural land decreased as the coverage of built-up areas increased. During the previous two decades, the world's population has grown exponentially, leading to an expansion of urban areas and decreased agricultural lands [67]. Several studies, including those by [68][69][70][71], found that as a city grew, the price of agricultural land climbed. The pressure on agricultural land has multiple dimensions. There was a transition from agricultural land to horticultural land in several areas. Additionally, agricultural land was transformed into a built-up area. Similarly, the forest area decreased steadily over the study period. As a result of the shift of thick forests into sparse forests and snow and glacier into barren and scrub regions, there was a favorable change in certain classes.
The predominant business of north Kashmir is horticulture, which uses most of the available land. Horticulture plantations have sprung up due to the incredible economic benefits of horticultural crops. According to [67,72,73], the primary causes of climate change are an increase in temperature and a decrease in precipitation, especially during peak cropping seasons. Increased evapotranspiration and the subsequent moisture shortage for crop production are the outcomes of these alterations, causing communities to change their farming practices. Shrinkage in forest cover also contributes to the expansion of scrubland because it leaves more area vulnerable to erosion and runoff. There was a similar decline in forest cover across the study area during the past two decades. Scrubby, dense, and scarce were the forests in the study area. Once thick with forest, the region now has only a few trees here and there. During the research, the area beneath water bodies has also decreased due to population growth and the expansion of built-up areas in the region. There has been a steady decline in the area covered by snow and glaciers as time goes on. The melting of glaciers is primarily due to a rise in temperature over the previous several decades, which has increased in areas covered by barren lands and scrub areas. In addition, the land-use shift in North-Western Himalayas, in general, and north Kashmir, in particular, has been susceptible to several interconnected variables. Temperature and precipitation extremes were significant drivers of land-use change in the study region. Due to the irregular precipitation patterns, farmers were obliged to adjust their land-use practices, which influenced agricultural output. Additionally, a rise in population and increased economic output from specific sectors helped maintain land-use change in the study area. Many studies have revealed that changes in land usage in the area have also been influenced by psychosocial characteristics [74]. In time, less and less of the land that had been left fallow for pasture was used for that purpose. According to [75], the increase in population has made it more challenging to keep fallow and pasture land in good condition. The rate of growth in the region's progress was meteoric.

Land Cover Prediction Using Transition Probabilities and Transition Matrix
The simulated and current LC dynamics indicated the North-Western Himalayas as a highly probable area of LC change, with horticulture and urban LC taking over other LC classes. The human population explosion continues to dominate urban LC at the cost of fallow and agricultural LC. The results of LC change simulation using the CA-Markov model shows that due to diverse climate conditions and higher economic returns, horticultural LC increased from 3.92% in 2000 to 8.43% in 2020 and is projected to reach 9.61% in 2030. Comparing the present, past, and future scenarios, horticulture, and urban areas have dominated the LC change. The changed areas in the spatial distribution show big urban hubs and horticultural development. Based on the CA-Markov model, these altered regions appear more prevalent in the north and west. The CA-Markov model appeared to provide logically credible findings for recognizing patterns. Only the models that successfully passed validation and calibration were used to build the future scenario. Similarly, the validity and certification of historical data satisfied the condition of an extended period for calibration. Furthermore, because LC is too dynamic and complex, spatial validating and calibrating of all changes across time and space can only be conducted using the CA-Markov model.

Conclusions
Since 2000, the North-Western Himalayan LULC has changed dramatically due to economic development and the influence of human activities. LC maps for 2000, 2010, and 2020 were created using Landsat 5 TM and Landsat 8 OLI images. The research area's LC structure was then simulated and projected using the CA-Markov model. The NWH have a high forest coverage, with areas of 1823, 1708, and 1695 km 2 in 2000, 2010, and 2020, respectively. Every year, construction land increased from 2000 to 2020. Changes in water, fallow land, and farming areas were linked to human activity. Land-use changes in the NWH from 2000 to 2020 were visible due to human activity. From 2000 to 2010, sand mining equipment was installed in the NWH, and a large amount of silt was deposited on the riverbank, showing a dramatic reduction in the water area. Although the modified region is enormous, and the percentage is modest, the wooded area is considerable. Timber harvesting and urban growth is the primary cause of change in the wooded area. The findings revealed that until 2030, the amount of woodland in the NWH will decline dramatically. We should pay attention to the quantity of woodland harvest in planning, considering the ecological functions of the forest. A few instances in the experimental method impacted the prediction outcomes. First, there were particular data collection challenges due to the research area's location.
Furthermore, the research region is in a hilly place where the ground varies in height, which influences the pixel value of the images and may lead to inaccurate classification results. Second, the setting of the suitability data-set had some effect on the LULC forecasts. Finally, there was some human influence on LULC types, particularly changes in building land. The CA-Markov model is a dynamic model that gives predictions in time and space and seems very practical. As a result, of the CA-Markov model simulations, the accuracy of the projected LC scenarios may be improved by enhancing the quality of the input data and changing associated parameters. Data Availability Statement: Your recommended statement did not apply on our article. Data supported this study is already mentioned in this manuscript.