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

Abstract

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 significant 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 classification of the study area using the maximum likelihood supervised classification. 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 first used to predict the land cover for 2020, which was then verified by the actual land cover of 2020 (Kappa coefficient 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 significant drivers of LULC changes were economic compulsions, climate variability, and increased human population. The analysis finding of the study highlighted that technical, financial, policy, or legislative initiatives are required to restore fragile NWH regions experiencing comparable consequences.

1. Introduction

In the past several decades, the population explosion has significantly altered the face of the globe, resulting in a massive change in land use and land cover (LULC). Changes in land use and land cover are considered one of the planetary boundaries that must be studied to maintain ecological balance and human existence [1,2,3,4,5]. The scientific study of LULC started in the late 1990s (Global Environment Change and Human Dimension Programme), giving rise to land use and land cover [6,7,8]. The terms “land use (LU)” and “land cover (LC) [9]” are used synonymously, but they represent two different schools of thought. LU change may be defined as physical, chemical, and biological change attributes that encompass the conversion of grazing, forest, and uncultivated land, to cropping land with alteration in fertilizer use, reclamation of degraded and polluted land, drainage improvement, installation of irrigation facilities, plantation, building farm dams, cultivation of exotic species [10,11,12,13,14]. LC denotes biological and physical cover on the land surface and includes water, bare soil, vegetation, and artificial constructions [1,15,16]. The LC has become a global concern in the present era [17,18]. The driving force of LC dynamics in space and time contributes to good infrastructure for natural resource sustainability, designing policies, and decision-making. On the other hand, land use is a predetermined avocation of land management planning used by land managers or other human beings to escapade the LC reflecting human activities such as crops, industries, forest plantations, residential buildings, airports, highways, logging, grazing, and many others [19,20,21,22,23,24,25,26,27,28]. The rapid and undesirable LULC changes and natural and anthropogenic activities have resulted in greenhouse gas emissions, habitat loss, biodiversity, desertification, urbanization, and wildfires, among others, affecting ecosystem processes and human existence [1,2,6,10,28]. These changes occur through the transgression of rural, urban, and forest lands, leading to severe environmental problems.

Policymakers, researchers, and scientists have been actively involved in determining the LULC change following several methods. Earlier, traditional procedures like cartography were followed. Conventional methods for LULC analysis were inadequate for multi-dimensional environmental studies because of the greater complexity of data and less accuracy of methods [29,30,31,32,33]. In this context, spatial data is paramount in determining quantitative and qualitative changes [34,35]. Over the past few years, advancements in remote sensing and GIS have enabled the repeated observation of the earth’s surface. With the enhanced capability of the sensors in terms of temporal frequency, spectral variability, and spatial resolution, small changes can be detected on the surface of the earth [9,36,37,38,39]. For larger areas, remote sensing and GIS are considered dynamic tools to obtain timely and meticulous instructions on the spatial distribution of LC changes. Numerous methods have been developed for change detection in land use and land cover. These include post-classification comparison, image regression, image change detection, image rationing, image ratio, conventional image differentiation, memory insertion, multi-data composite, multivariate image classification, vegetative difference indexing, manual on-screen digitization component, and land use and land cover modeling [3,40,41]. The commonly used models for observing and recording land use and land cover changes are statistical, analytical equation-based, cellular, evolutionary, hybrid, Markow, multi-agent, and expert system models [42,43,44]. The most commonly used land use and land cover models are cellular and agent-based.

The studies on the LULC are of paramount importance for understanding the factors responsible for these changes, including climate change, economic viability, ecosystem change, and the growth of urban areas. The harmful impact of the unplanned LULC can be reduced through suitable land management strategies, which need knowledge about the present, past, and future scenarios of LULC. Therefore, there arises a need for time series data for LULC. This research was carried out in the Baramulla district of Jammu and Kashmir’s Union Territory, which has been experiencing unchecked urban growth since 2000. This study examines how the distinct LULC processes can be accurately simulated and predicted by factoring in their underlying geographical determinants. To 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.

2. Material and Methods

2.1. 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²) 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. Site characters and natural vegetation of the study area.

Land Cover	Latitude and Longitude	Elevation (m.a.s.l.)	Topography	Slope (%)	Depth of Soil (cm)	Natural Vegetation
Agriculture	34°10′19″ N 74°31′02″ E	1983	Undulating	3–8	0–179	Pinus spp., Ulmus spp., Populus spp., Salix spp., Fir spp., Berberis spp., Aaicheria spp.
Horticulture	34° 12′57″ N 74° 21′49″ E	2385	Rolling	8–16	0–83	Pinus spp., Ulmus spp., Populus spp., Wild grass spp., Walnut spp., Celtis spp., Aaicheria spp.
Agro-Forestry	34°15′50″ N 74°18′18″ E	2162	Foot hills	16–25	0–188	Pinus spp., Ciderus spp., Populus spp., Ailanthus spp., Walnut spp., Urtica spp., Aaicheria spp., Rumex spp.
Fallow Land	34°2′32″ N 74°14′06″ E	2110	Rolling	8–16	0–114	Populus spp., Salix spp, Walnut spp, Taraxicum spp., Malwa spp., Berberis spp., Cotoneaster spp., Aliesthus spp.

).

2.2. 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:

$$\Delta \mathrm{A} (\%)=\frac{{At}_2-{At}_1}{{At}_1}$$

(1)

A (%) = percentage of difference in land use and land cover type between the initial and final value, At₁ is the initial area of land use and land cover, and At₂ 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:

$$\frac{\Delta \mathrm{R}\left(\mathrm{ha}\right)}{\mathrm{Yr}}=\frac{\mathrm{C}-\mathrm{P}}{\mathrm{T}}$$

(2)

where;

ΔR = rate of change, C = current area coverage of LC, P = previous area coverage of LC, T = time interval between C and P

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. Area proportion of LC types in 2000, 2010, and 2020.

Land Cover	2000		2010		2020
	km2	%	km2	%	km2	%
Water	74.00	1.06	68.00	0.98	54.00	0.77
Snow	400.00	5.74	357.00	5.12	312.00	4.48
Vegetation	1280.00	18.36	1315.00	18.86	1365.00	19.58
Forest	1823.00	26.15	1708.00	24.50	1695.00	24.31
Agriculture	1555.00	22.30	1484.00	21.29	1397.00	20.04
Horticulture	273.00	3.92	402.00	5.77	588.00	8.43
Urban	677.00	9.71	918.00	13.17	984.00	14.11
Fallow land	890.00	12.77	720.00	10.33	577.00	8.28
Total	6972	100	6972	100	6972	100

). 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].

2.3. 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.

2.3.1. 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.

$${\mathrm{L}}_{c+1}=f({\mathrm{L}}_c+N)$$

(3)

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.

2.3.2. 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].

$$=\left[\begin{array}{ccc}{\mathrm{L}}_{11}& \dots & {\mathrm{L}}_{1\mathrm{n}}\\ ⋮& ⋮& ⋮\\ {\mathrm{L}}_{\mathrm{n}1}& \dots & {\mathrm{L}}_{\mathrm{mn}}\end{array}\right]$$

(4)

$$\left(0\le {\mathrm{P}}_{\mathrm{ij}}<1 \mathrm{and}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{P}}_{\mathrm{ij}}=1, \mathrm{i},\mathrm{j}=1, 2, \dots , \mathrm{n}\right)$$

(5)

$${\mathrm{L}}_{\left(c+1\right)}={\mathrm{P}}_{\mathrm{ij}}\times {\mathrm{L}}_c$$

(6)

where n is the LC type number, L is LC type.

2.3.3. 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 decision-making 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.

2.3.4. 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 information (Kno) [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.

2.3.5. 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].

$${\mathrm{L}}_{\mathrm{i}}=\frac{{\mathrm{LA}}_{\left(\mathrm{i},\mathrm{t}1\right)}-{\mathrm{ULA}}_{\mathrm{i}}}{{\mathrm{LA}}_{\left(\mathrm{i},\mathrm{t}1\right)}}\times \frac{1}{{\mathrm{c}}_2-{\mathrm{c}}_1}\times 100\%$$

(7)

$$\mathrm{L}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left\{{\mathrm{LA}}_{\left(\mathrm{i},\mathrm{t}1\right)}-{\mathrm{ULA}}_{\mathrm{i}}\right\}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{LA}}_{\left(\mathrm{i},\mathrm{t}1\right)}}\times \frac{1}{{\mathrm{c}}_2-{\mathrm{c}}_1}\times 100\%$$

(8)

$$\mathrm{DOR}=\mathrm{TOR}+\mathrm{TIR} 9$$

(9)

$$\mathrm{TOR}=\frac{{\mathrm{AP}}_{\left(\mathrm{i},\mathrm{t}1\right)}-{\mathrm{ULA}}_{\mathrm{i}}}{{\mathrm{LA}}_{\left(\mathrm{i},\mathrm{t}1\right)}}\times \frac{1}{{\mathrm{c}}_2-{\mathrm{c}}_1}\times 100\%$$

(10)

$$\mathrm{TIR}=\frac{{\mathrm{AP}}_{\left(\mathrm{i},\mathrm{t}2\right)}-{\mathrm{ULA}}_{\mathrm{i}}}{{\mathrm{LA}}_{\left(\mathrm{i},\mathrm{t}1\right)}}\times \frac{1}{{\mathrm{c}}_2-{\mathrm{c}}_1}\times 100\%$$

(11)

where L_i and L are the dynamic degrees of single and comprehensive LULC dynamic degrees, LA (i, c₁) is the land-use area of a particular LC type at an initial stage, and LA (I, c₂) is the land-use area of a particular LULC type at a later stage, ULA_i is the area not changed, c₁ and c₂ 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.

3. Results

3.1. 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. Accuracy assessment reports of the image classification at the different temporal regimes.

LC	2000		2010		2020
	Producer’s Accuracy	User’s Accuracy	Producer’s Accuracy	User’s Accuracy	Producer’s Accuracy	User’s Accuracy
Water	76.19	81	79.24	86.14	84.1	89.32
Snow	72.9	75.2	76.31	78.23	81.3	89.6
Vegetation	88.56	82.25	85.14	89.35	82.4	77.63
Forest	94.03	97.44	89.26	93.54	91.54	96.37
Agriculture	81.43	83.89	92.45	97.36	89.56	91.18
Horticulture	86.24	89.24	93.75	98.82	95.52	98.91
Urban	74.35	77.38	86.48	89.93	83.09	88.62
Fallow land	80.68	78.57	84.9	89.64	94.34	95.35
Kappa coefficient	81.79	83.12	85.94	90.37	87.73	90.87
Overall accuracy	81.13		78.09		81.0

).

3.2. Land Cover Transformation Dynamics

Table 4. Area changes in each LC type from 2000 to 2010, 2010 to 2020, and 2000 to 2020.

LC	2000 to 2010		2010 to 2020		2000 to 2020
	km2	%	km2	%	km2	%
Water	6.00	8.11	14.00	20.59	20.00	37.04
Snow	43.00	10.75	45.00	12.61	88.00	28.21
Vegetation	−35.00	−2.73	−50.00	−3.80	−85.00	−6.23
Forest	−115.00	−6.31	−13.00	−0.76	−128.00	−7.55
Agriculture	71.00	4.57	87.00	5.86	158.00	11.31
Horticulture	129.00	47.25	186.00	46.27	315.00	53.57
Urban	241.00	35.60	66.00	7.19	307.00	−1.20
Fallow land	170.00	19.10	143.00	19.86	313.00	54.25

depicts the change in LULC during 2000–2010, 2010–2020, and 2000–2020 with specific gain or loss in each LULC type. The statistical data in

Table 5. The change rate of LC during 2000–2010.

LC Type	Total (Area/km2)	Unchanged (Area/km2)	Transfer Area (Area/km2)	Transfer Rate (Area/km2)	Gain Area (Area/km2)	Gain Rate (Area/km2)	Rate of Change (Area/km2)	Dynamic Degree %
Water	74.00	58.19	15.81	2.14	9.81	1.33	3.46	3.53
Snow	400.00	330.47	69.53	1.74	26.53	0.66	2.40	0.91
Vegetation	1280.00	1235.56	44.44	0.35	79.44	0.62	0.97	0.12
Forest	1823.00	1678.09	144.91	0.79	29.91	0.16	0.96	0.32
Agriculture	1555.00	1467.44	87.56	0.56	16.56	0.11	0.67	0.40
Horticulture	273.00	250.09	22.91	0.84	151.91	5.56	6.40	0.42
Urban	677.00	602.98	74.02	1.09	315.02	4.65	5.75	0.18
Fallow land	890.00	654.72	235.28	2.64	65.28	0.73	3.38	0.52

and

Table 6. The rate of LULC change during the 2010– 2020 period.

LC	Total (Area/km2)	Unchanged (Area/km2)	Transfer Area (Area/km2)	Transfer Rate (Area/km2)	Gain Area (Area/km2)	Gain Rate (Area/km2)	Rate of Change (Area/km2)	Dynamic Degree (%)
Water	68.00	48.79	19.21	2.83	5.21	0.77	3.59	6.89
Snow	357.00	286.81	70.19	1.97	25.19	0.71	2.67	1.06
Vegetation	1315.00	1299.54	15.46	0.12	65.46	0.50	0.62	0.09
Forest	1708.00	1662.97	45.03	0.26	32.03	0.19	0.45	0.14
Agriculture	1484.00	1313.44	170.56	1.15	83.56	0.56	1.71	0.20
Horticulture	402.00	370.1	31.90	0.79	217.90	5.42	6.21	0.29
Urban	918.00	734.34	183.66	2.00	249.66	2.72	4.72	0.19
Fallow land	720.00	515	205.00	2.85	62.00	0.86	3.71	0.60

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². The amount of horticultural land-use area converted from another LC was 129 km². 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², accounting for 6.31%. From 2010 to 2020, the horticultural land cover had the greatest change, with an area of 186 km² 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², 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². 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² and 307 km², respectively. Horticultural land was transformed from fallow land, followed by urban, water, snow, forest, agriculture, and vegetation.

Table 5 and Table 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⁻¹ 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⁻¹ 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⁻¹. The lowest transformation rate was observed in vegetation with a value of 0.35 and 0.12 kma⁻¹. 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⁻¹. The lowest gain rate occurred in a forest area with a value of 0.6 and 0.19 kma⁻¹. 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.

3.3. Area Transfer and Probability Matrix of LULC

Each row in

Table 7. Area transfer matrix from 2000–2010 (km²).

LC	Water	Snow	Vegetation	Forest	Agriculture	Horticulture	Urban	Fallow Land	Total
Water	58.19	0.18	7.76	1.23	2.19	0.34	3.00	1.11	74.00
Snow	15.67	330.47	8.97	12.56	0.15	0.89	0.00	31.29	400.00
Vegetation	2.12	0.45	1235.56	8.89	3.45	12.78	14.13	2.62	1280.00
Forest	2.68	15.98	20.09	1678.09	30.88	57.09	16.09	2.10	1823.00
Agriculture	2.09	0.34	1.03	0.09	1467.44	69.00	13.56	1.45	1555.00
Horticulture	0.05	0.11	1.56	1.40	3.12	250.04	14.99	1.73	273.00
Urban	0.55	0.13	4.01	5.17	18.43	42.28	602.98	3.45	677.00
Fallow land	1.23	0.76	3.46	38.90	50.77	70.87	69.30	654.71	890.00
Total	82.58	348.42	1282.44	1746.33	1576.43	503.29	734.05	698.46	6972.00

and

Table 8. Transfer Probability matrix during 2000–2010 (km²).

LC	Water	Snow	Vegetation	Forest	Agriculture	Horticulture	Urban	Fallow Land
Water	59.70	3.14	4.05	1.89	8.12	7.60	5.40	10.10
Snow	10.03	48.6	5.56	14.92	1.90	0.94	3.00	15.05
Vegetation	2.60	1.45	38.79	15.80	6.08	20.09	10.06	5.13
Forest	0.31	1.8	4.89	65.48	3.13	13.65	6.13	4.61
Agriculture	0.45	0.21	2.34	4.56	51.80	24.60	11.90	4.14
Horticulture	0.11	0.07	1.78	2.54	1.67	89.09	3.79	0.95
Urban	0.06	0.34	2.56	3.45	4.13	9.79	78.98	0.69
Fallow land	0.67	0.35	4.89	7.90	10.62	16.57	24.89	34.11

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², 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², 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², and the snow area decreased by 43 km², accounting for 8.11 and 10.75%. The water and snow area are transformed into vegetation and forest LC.

Table 9. Area transfer matrix from 2010–2020 (km²).

	Water	Snow	Vegetation	Forest	Agriculture	Horticulture	Urban	Fallow Land	Total
Water	48.79	0.24	5.45	2.78	3.12	0.67	1.45	5.50	68.00
Snow	11.11	286.81	14.45	17.89	0.87	1.22	0.08	24.57	357.00
Vegetation	0.56	0.09	1299.54	4.67	1.23	5.67	3.34	0.46	1315.00
Forest	0.12	1.32	2.34	1662.97	12.78	21.09	5.67	1.83	1708.00
Agriculture	0.23	0.067	8.98	12.34	1313.08	98.99	37.98	12.56	1484.00
Horticulture	0.56	0.012	1.56	4.67	2.30	370.1	21.89	1.47	402.00
Urban	1.14	0.01	12.30	6.23	14.13	36.09	846.34	2.90	918.00
Fallow land	3.89	5.3	12.34	20.09	39.36	93.21	34.70	515	720.00
Total	66.40	293.85	1356.96	1731.64	1386.87	627.04	951.45	564.29	6972

and

Table 10. Transfer Probability matrix during 2010–2020 (km²).

	Water	Snow	Vegetation	Forest	Agriculture	Horticulture	Urban	Fallow Land
Water	53.70	2.43	7.89	2.09	9.18	2.21	12.40	10.10
Snow	9.89	34.6	10.01	10.09	3.79	0.24	1.33	30.05
Vegetation	3.17	2.31	49.89	10.23	6.24	17.09	7.06	4.01
Forest	1.43	2.19	5.80	60.16	7.88	17.21	3.12	2.21
Agriculture	1.80	0.08	0.52	2.12	42.69	36.60	15.05	1.14
Horticulture	0.17	0.01	2.48	1.09	3.67	85.09	7.24	0.25
Urban	1.10	0.06	1.50	4.67	3.45	13.45	74.97	0.80
Fallow land	1.63	0.57	1.22	3.12	8.00	25.57	17.37	42.52

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², 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², which accounts for 13.77%, and its expansion is due to agriculture horticulture land transfer-in, agriculture land transfer-in, vegetation land transfer-in, fallow forest land transfer-in, and fallow land transfer-in with transferring area of 36.09, 14.13, 12.30, 6.23, and 2.90%.

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. Projection of LULC of the area in 2030.

Land Cover Class	Projected Area Coverage (km2)	Area Proportion (%)
Water	63	0.90
Snow	317	4.55
Vegetation	1345	19.29
Forest	1630	23.38
Agriculture	1393	19.98
Horticulture	670	9.61
Urban	1034	14.83
Fallow land	520	7.46
Total	6972	100

).

4. Discussion

4.1. 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 Table 3 and Table 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.

4.2. 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.

5. 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² 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.