Land Use and Land Cover Change Modeling and Future Potential Landscape Risk Assessment Using Markov-CA Model and Analytical Hierarchy Process

Abstract

Land use and land cover change (LULCC) has directly played an important role in the observed climate change. In this paper, we considered Dujiangyan City and its environs (DCEN) to study the future scenario in the years 2025, 2030, and 2040 based on the 2018 simulation results from 2007 and 2018 LULC maps. This study evaluates the spatial and temporal variations of future LULCC, including the future potential landscape risk (FPLR) area of the 2008 great (8.0 M_w) earthquake of south-west China. The Cellular automata–Markov chain (CA-Markov) model and multicriteria based analytical hierarchy process (MC-AHP) approach have been considered using the integration of remote sensing and GIS techniques. The analysis shows future LULC scenario in the years 2025, 2030, and 2040 along with the FPLR pattern. Based on the results of the future LULCC and FPLR scenarios, we have provided suggestions for the development in the close proximity of the fault lines for the future strong magnitude earthquakes. Our results suggest a better and safe planning approach in the Belt and Road Corridor (BRC) of China to control future Silk-Road Disaster, which will also be useful to urban planners for urban development in a safe and sustainable manner.

1. Introduction

The land use land cover change and its modeling (LULCC-M) approach has recently been considered by the scientific community to observe environmental changes at level 4, such as local, regional, national, and global. To achieve the Sustainable Development Goals framework designed by the United Nations (UN) 2030, Sustainable Development Goals (SDGs) have been considered worldwide. Changes above level 4 are heavily impacted by natural hazards, i.e., earthquakes, floods, landslides, and forest fires. The term ‘land use and land cover changes’ (LULCC) refers to human modifications of the terrestrial surface of the earth, as well as the study of land surface change [1], and the term ‘transition’ defines the process of changing, or the change of something from one form or state to another [2]. The study of LULC change is considered as an important aspect to observe parameters which are responsible for overall changes [3,4,5,6,7,8,9,10,11,12,13,14,15,16], landscape changes [17,18,19,20,21,22], landscape fragmentation [23,24], changes in ecosystems [25,26], climate changes [27], growth of urban areas [28,29], sustainable development [30,31], and environment and risk evaluation of earthquake-affected cities [16]. Sui and Ming’s [31] idea of LULCC on sustainable development in China has attracted the attention of scientists and planners worldwide. However, LULCC at spatiotemporal scales in metropolitan areas has been carried out worldwide in earthquake-affected cities, watersheds, coastal areas, and wildlife sanctuary areas using multitemporal remote sensing satellite data [16,32,33,34,35,36,37,38,39,40,41,42,43].

To support LULCC, land-use change models (LUCM) are considered as effective tools to analyze the causes and consequences of land-use change. Moreover, the LULC model is found to be effective when predicting the future state and spatial distribution of LULC using the gained knowledge from previous years [44]. Many scientists, such as Singh et al. [45,46,47,48,49], Dey and Singh [50], Okada et al. [51], and Dey et al. [52], have conducted studies related to surface, ocean, and atmospheric parameters associated with the Gujarat earthquake of 26 January 2001. Singh and Singh [53] also performed lineament mapping associated with the changes in stress in the epicentral areas of the Gujarat earthquake using remote sensing data. In addition, group of scientists carried out a detailed analysis of damages and numerous parameter changes associated with the deadly 2008 Wenchuan earthquake which occurred in south-west (SW) China [54,55,56,57,58,59,60,61]. The remote sensing and GIS techniques have been widely used to assess LULCC, especially in inaccessible mountainous regions [62,63,64,65]. Moreover, high-resolution remote sensing images have also been used in efficient image segmentation [66] and knowledge transfer to automatically update LULCC databases [67] in different areas of China. In connection with the LULCC study, research on the lineament structural change caused by earthquakes in South America and China has been carried out using ASTER (Terra) [68] and Landsat (5-TM and 8-OLI) [16] satellite data, respectively.

Numerous methods, such as mathematical-equation-based, spatiotemporal modeling [69], system dynamic simulation [70], statistical, cellular and hybrid models [71], cellular and agent-based models or a hybrid of the two [72], and the cellular automata–Markov chain (CA-Markov) model [44,73,74,75,76], have been utilized in different research. CA-Markov initial conditions, parameterization of the model, calculations of the transition probabilities, and determination of the neighborhood rules were defined by the remote sensing and GIS datasets [77,78,79,80,81]. The CA-Markov model is one of the most ideal and widely accepted methods for LULCC modeling because it considers ‘t-1’ to ‘t’ to project probabilities of LULCC for the future date ‘t+1’ [44,82,83]. The probabilities are generated based on past and future changes [44]. The CA-Markov model has the ability to simulate changes in different LULC and can possibly simulate the transition from one category of LULC change to another [44]. However, a combined CA-Markov model to simulate future LULCC by integrating natural and socioeconomic data is still a challenge due to the different datasets.

In Sichuan province, many earth observation data have been used for land management, monitoring, and degradation studies [84,85,86,87], especially after the Wenchuan earthquake. Dujiangyan City and its environs (DCEN) are located at the southeast of the Tibet plateau under Sichuan province. Asia’s Water Tower (AWT) suffers from frequent geological disasters, and therefore, the region needs close attention. In the past, DCEN were not studied in detail, but by observing the area in 2018, Nath et al. [16], studied major landscape changes in the area first time using Landsat datasets. The area is close to the epicenter of the deadly Wenchuan earthquake of 12 May 2008, and they observed the LULCC, including the landscape risk patterns (LRP) that developed in the 2007–2018 time period. However, in the DCEN region, risk evaluation studies are limited. The present study focuses in LULCC for the period 2007–2018 which will help to make plans and take steps to improve ecological environments, which are a threat due to future changes in LULC. Therefore, the present study is an important effort to closely monitor the DCEN area so that appropriate action can be considered. In connection with this, future land use and land cover change (FLULCC) and future potential landscape risk (FPLR) patterns for the years 2025, 2030, and 2040 were observed based on the cellular automata–Markov chain (CA-Markov) modeling and multicriteria-based analytical hierarchy process (MC-AHP) approach. LULCC modeling and FPLR assessment was carried out by integrating remote sensing and GIS techniques developed for the first time in the DCEN area. However, this study considers it a unique approach to model the DCEN area closely, focusing on FLULCC along with the FPLR scenarios up to 2040.

Comparison and analysis of LULC in the distant past and in recent times is important in order to understand the changes in spatiotemporal LULC [88,89,90,91,92]. After the earthquake, the China government began to rebuild the city in 2008–2010, where the survivors along with new residents from outside areas started living all together without having any prior or aftermath knowledge, including any further presence or not of risk in their locality. In this research, two multitemporal Landsat-image-based LULC maps of 2007 and 2018 (adopted with permission from Nath et al. [16]) were considered for the simulation study of 2018. The year 2007 was considered for actual existing LULC prior to the earthquake (BEQ) in 2008 (8.0 M_w), and January 2018 is considered for recent LULC changes and the LRP of the study area. The disaster chain effects are the main properties of knowing the landscape risk. Based on the method followed by Nath et al. [16], in this study, FPLR scenarios were also developed using Remote sensing and GIS tools. Moreover, multiple parameters such as earthquake epicenter locations, local geology, fault lines, geomorphological features, and Lineament density (LD) of 2018 and future projected LULC datasets of 2025, 2030, and 2040 were considered to know the FPLR scenarios in the study area, where few parameters are likely to change with time.

The present study focuses on quantitative analysis of spatiotemporal changes of future LULC changes in 2025, 2030, and 2040 based on the 2018 LULC simulation and further considers recent LR 2018 data with ancillary layers to identify FPLR areas across the DCEN. In this study, examples are shared for each future LULC map preparation, referring to 2018 as the earlier example and 2025 as the recent one, and so on. However, the broad objectives of this study were to examine three issues: (1) quantitatively evaluating FLULC changes in the DCEN in considering 2018 simulation data, (2) quantitatively evaluating future LULCC in different time domain matrices’ data such as comparing earlier LULC datasets with the recent one (2018–2025, 2025–2030, 2030–2040, and 2018–2040) used for the DCEN regions; and (3) identifying and evaluating FPLR areas in consideration of important multiple environmental parameters. This study will be useful in understanding the above three objectives for the landscape planning and environmental management in this area in an effective and sustainable manner. Finally, it will be contributing to achieving the Goal 11 (“Sustainable Cities and Communites”) of the United Nations (UN) 2030 Sustainable Development Goals (SDGs) to make cities inclusive, safe, resilient and sustainable [93].

2. Materials and Methods

2.1. Study Area

DCEN is a county-level city located in the Sichuan province of China, which covers an area of 318.30 km². The city, which suffered a deadly earthquake in 2008, is now a UNESCO world heritage site and underwent massive renovations in 2013. It shares a border with southern Nagawa Tibet in the north, due to which it receives many tourists. The geography of the study site is mountainous with many forests, rivers, bridges, and urban infrastructures [16]. The study area of DCEN (30°53′0″ N to 31°4′30″ N and 103°27′5″ E to 103°43′0″ E) is an earthquake-affected city in SW China, part of the provincial capital of Chengdu of Sichuan province (Figure 1).

This study area of DCEN was first adopted by Nath et al. [16] and has been considered in the present study with prior permission to investigate the LULCC modeling and FPLR assessment. DCEN has a humid subtropical climate (Köppen Cwa) with cool, dry winters and hot, very wet summers. DCEN has 17 towns, including two townships. According to Google Earth, the elevation and slope in DCEN mostly varies from 668 to 1861 m from above mean sea level (AMSL), and the maximum slope and average slope are 67.9% and 18.9%, respectively. The study area has suffered floods and landslide disasters in the past, and greater earthquakes and aftershocks occurred in 2008–2009 and afterwards.

2.2. Data

The remote sensing data used in this study were based on the Landsat 5 thematic mapper (TM) and Landsat 8 operation land imager (OLI). The spectral resolution and bands of the data were used in this study (

Table 1. Details of data used in the present study.

SL. No.	Satellite Sensors	Year	Image Acquisition Time	Path/Row	Resolution (m)	Image Type
1	Landsat 5 TM	2007	18 September 2007	130/38	30	Level-1 GeoTIFF
2	Landsat 8 OLI	2018	19 January 2018	130/38	30	Level-1 GeoTIFF

Note: In the 2nd column TM refers to thematic mapper, and OLI refers to operation land imager. Source: USGS Earth Explorer Landsat archive (https://earthexplorer.usgs.gov) [94].

). Two Landsat satellite images, Landsat 5 TM in 2007 and Landsat 8 OLI in 2018, were considered in the present study. However, it is worth noting that in this study, the two datasets were originally obtained from the United States Geological Survey (USGS) Landsat archive: 2007 and 2018 [94]. The choice of those images was made to observe future LULCC modeling and FPLR assessment.

The image covers less than 5% cloud in original satellite images, 2% in 2007, 4.65% and 2018 (4.65%). Due to the presence of cloud cover in the epicentral area, we only considered the area near the 2008 earthquake epicenter (8.0 M_w). The area of interest (AOI) is demarcated in the satellite images with the help of the polygon feature using the ArcGIS 10.6 software (ESRI, Redlands, USA) environment. The two Landsat images were geometrically corrected and projected to a WGS 84 (World Geodetic System 1984) into a Universal Transverse Mercator (UTM) Zone 48N coordinate system. Detailed preprocessing of the Landsat 5 TM and 8 OLI images was carried out by Nath et al. [16] to get surface reflectance images. The Landsat classified images of LULC of 2007 and 2018 (adopted from Nath et al. [16]) were used in the present study with permission and with the color code changed.

In addition, a shuttle radar topographic mission (SRTM) with 30m resolution was used for the digital elevation model (DEM) (Figure 1), followed by a derived topographic feature (i.e., slope). Moreover, the frequency of M_w 4.0 and higher earthquakes in the epicenter and its vicinity regions is shown for the periods 12 May 2008 to 8 September 2019 in Figure 2 (modified from Nath et al. [16]).

All the maps have been prepared at a 1: 1,25,000 scale except Figure 2, which was prepared at a 1: 35,00,000 scale. A combination of geological maps, fault lines, epicenter points, FLULC of 2018, 2025, 2030, and 2040, and the lineament density (LD) map of 2018 were used to identify the FPLR areas.

2.3. Methodology

The integrated research methodological workflow requires four main steps (Figure 3). Step 1: Image preprocessing of two multitemporal remote sensing images, Landsat 5 TM and Landsat 8 OLI, image subsetting based on AOI and modification of multiple layers using a reclassifying operator, image classification using MLC techniques and accuracy assessment, and then model validation and simulated LULC mapping of 2018 generation with TerrSet 18.2 software based on the two images. Step 2: Creation of FLULC maps based on the CA-Markov model first, then transition matrix table creation based on different time periods with dynamic degree estimation, including overall change. Step 3: PCA image creation based on the Landsat 8 OLI image and the lineament-based lineament density (LD) map adopted from Nath et al. [16], multiple layer integrations, including individual FLULC maps and weighted-score-based FPLR risk maps, and thematic classification performed, as shown in Figure 3.

2.3.1. LULCC Modeling Using Land Change Modeler

In the present study, we used the Clark Lab’s land change modeler (LCM) of TerrSet 18.21 version software to develop FLULC maps based on the preprocessed classified historical satellite images of 2007 and 2018 (adopted from Nath et al. [16]). TerrSet is an integrated geospatial software system widely used by the scientific community to monitor and model the earth system processes for sustainable development. However, the LCM was especially developed for the vertical application to analyze LULCC and empirical modeling and project future changes. These two classified maps were used to model the LULC of the study area for the future time periods (i.e., 2025, 2030, and 2040). In addition, multiple GIS data layers have been considered and performed subsetting based on AOI along with multiple layers using the reclassifying operator of ArcGIS 10.6 version software.

In order to run the future predictions, the integrated CA-Markov approach [44] was used. Based on the 2007 and 2018 images, the transition probability matrix (TPM) and transition probability areas (TPA) were prepared from the Markov model. The transition suitability image (TSI) was also produced. Finally, the TPA, TPM, and TSI were integrated in the CA-Markov model. A 5*5 contiguity filter was applied with 5 iterations to model the LULC for the years 2025, 2030, and 2040.

In addition, the model further needed to be validated first before running the simulation. Since the actual LULC image of 2018 was already available, adopted from Nath et al. 2018 [16], and was compared with the simulated LULC map of 2018 to evaluate its final agreement, the model validation was carried out using the Markov chain TPM of LULC classes for the periods 2007–2018 (shown in

Table 2. Markov chain transition probability matrix (TPM) of land use land cover (LULC) classes in 2018 (simulated) by percentage based on 2007 and 2018 LULC data.

2018
	FLULC Classes	BU	F	AG	RA	WB	BL
	BU	0.8987	0.0025	0.0000	0.0330	0.0148	0.0510
	F	0.0000	0.9894	0.0000	0.0098	0.0003	0.0006
2007	AG	0.2292	0.0683	0.0000	0.1011	0.0000	0.6014
	RA	0.4744	0.0000	0.0000	0.2385	0.0000	0.2871
	WB	0.0000	0.0000	0.0000	0.0303	0.9697	0.0000
	BL	0.0000	0.0228	0.5693	0.0000	0.0000	0.4079

). The validation was done using the Kappa index of agreement and chi-square test statistics. The simulated 2018 image was used for the FLULCC prediction. Moreover, the TPM of FLULC classes in each time periods (such as, 2018–2025, 2025–2030, and 2030–2040) was obtained using the Markov model.

2.3.2. FLULC Dynamic Degree Estimation and Transition Matrices Computation Method

To estimate the dynamic degree (DD) changes between the various types of FLULC, i.e., to know whether there has been any loss in the area or gain in certain LULC types when comparing images for different periods, for a single FLULC type, we used the DD model to represent the spatiotemporal characteristics of FLULC change. The dynamic estimation can be calculated using the approach of Liu et al. [7,8,9]:

$$D=\frac{\frac{Ab-Aa}{Aa}}{T}\times 100\%$$

(1)

where D represents the DD model, which refers to rate of change; A_a is the area in the initial year; i is the area in the terminal year; and T is the temporal scale. In our case, the time comparisons are 7, 5, 10, and 22 years, respectively.

The DD model finally serves to generate statistical data in a table with DD (%) and Gain and Loss (%). Furthermore, transition matrix (TM) tables of 2018–2025, 2025–2030, 2030–2040, and 2018–2040 were computed based on the simulated LULC data of 2018, which were generated from the original 2007 and 2018 image results. Therefore, based on the 2018 simulated LULC results, similar statistical results for the different time periods of 2018–2025, 2025–2030, 2030–2040, and 2018–2040 are prepared in this study. These data (vector format) were analyzed and conversion was carried out using ArcGIS 10.6 software, which allows the identification of the areas where changes will occur, and those LULC were not altered in the FLULC. Furthermore, all of the results obtained using GIS software were exported to a text file and the databases were later used for statistical analysis.

2.3.3. FPLR Evaluation Method

In this study, the LR evaluation methodological concept adopted from Nath et al. [16], Richards and Jia [92], Anderson et al. [95], and Singh et al. [96] was used to determine the FPLR exposure based on the 2018 data and multiple features of the earth surface, including the projected LULC data of the years 2025, 2030, and 2040. This type of integrated assessment results will help planners to know about the potential threats in the coming decades in their areas, which might affected people, property, and the environment. Geology is one of the important parameters used for earthquake risk evaluation in DCEN for urban development. The concept of Tudes and Yigiter [97] which says that “development will take place where there is low risk of earthquake” was utilized for the suitability analysis of the city of Adana in Turkey. Thereafter, by modifying this concept, Nath et al. [16] focused on the 2008 earthquake which affected Dujiangyan City, where the LR areas with defined categories were identified and an observation of cities’ growth towards an active fault line was suggested to control further development. The concept of LR mapping and procedural steps was already discussed by Nath et al. [16]; therefore, we have not included such a discussion here. The lineament map from the satellite image was prepared using PCI image, followed by the corresponding lineament density (LD) map of 2018, which was further used in the FPLR assessment, subject to change with time.

Recently, the multicriteria evaluation (MCE) method has been used to assess and aggregate weighted maps (criterion) based on expert knowledge of factor influence and interactions with LULC [92,93]. Each FPLR map is derived using MCE by combining the information from multiple criteria to form a single index of evaluation [92,96,97].

Multiple parameters, such as the future LULC maps of 2025, 2030, and 2040, the 2018 lineament feature, the 2018 LD map, fault lines, and the epicenter, overlaid along with the local geology map, slope, aspect, and DEM, are used as reference datasets in integrated mapping. Moreover, the geology map was adopted from Nath et al. [16], which was modified and re-attributed by them as zone I–IV according to its existing category [98]. In this study, based on the above important layers, integrated overlay maps were prepared in three different time frames, 2025, 2030, and 2040, by following the previous method [16]. In the next stage, multiple buffer creations (i.e., 5, 10, 15, 20, and 25 km distance) and weightage (5, 4, 3, 2, 1) were given for each layer, such as important places, past earthquake epicenters, and fault lines, and along with LD-2018 and LULC-2025, 2030, 2040 data sets were considered further for FPLR evaluations. The weighted score value assigned was similar to the score followed by Nath et al. [16], i.e., “higher the risk value” when “buffer distance is lower”, and “lower the risk value” when “buffer distance is higher”. FLULC class was designated with risk priority from very high to very low followed by a previously adopted weighted score (Nath et al. 2018), such as, BU=1, RA=1, F=2, WB=3, BL=4, and AG=5. In this aspect, FLULC classes were designated as BU—built-up area; F—forest area; AG—agricultural area; RA—reconstructed area; WB—water bodies; and BL—barren land. This order was finally adopted and followed in the FLULC, FLULCC mapping, and FLULCC matrix table generation.

Details of the meaning and individual weighted assigned score of each layer are given in Supplementary Table S1a–g. This table is based on certain conditions and should be interpreted as: If the category of FPLR is ‘very high risk’ and ‘high risk’, it is close to the earthquake epicenter. On the other hand, ‘medium risk’ represents close to the high-risk area. Similarly, when the value is low, it indicates ‘low risk’ to ‘very low risk’, i.e., the risk is low because the area is far away from the past epicenters and geological fault lines. Reclassification was performed for all the layers with the aid of the field calculator of ArcGIS 10.6 software. In the next section, a new total risk weighted value field is created and weighted score is assigned automatically in the integrated database table for the FPLR maps of 2025, 2030, and 2040 in a consecutive manner. The geoprocessing wizard of the ‘union’ operator function of ArcGIS 10.6 has been applied to combine the data values of all layers. The total FPLR was defined by the integration of all weighted values as shown in Equation (2).

$$\begin{array}{ll}\mathrm{FPLR}\left[\mathrm{T}\_\mathrm{Weight}\right]& = \left[\mathrm{IP}\_\mathrm{B}\_\mathrm{Weight}\right]+\left[\mathrm{EQ}\_\mathrm{B}\_\mathrm{Weight}\right]+\left[\mathrm{F}\_\mathrm{B}\_\mathrm{Weight}\right]+\left[\mathrm{Risk}\_\mathrm{Geol}\right]\\ & + [\mathrm{LD}2018\_\mathrm{Risk}]+[\mathrm{FLULC}\left(\mathrm{i}.\mathrm{e}.2020,2030,\mathrm{and}2040\_\mathrm{Risk}\right]\end{array}$$

(2)

where FPLR[T_Weight] refers to future potential landscape risk (total weight); [IP_B_Weight] refers to important places buffer weight; [EQ_B_Weight] refers to earthquake point buffer weight; [F_B_Weight] refers to fault-based buffer weight; [Risk_Geol] refers to risk factor of geological zone; [LD2018_Risk] refers to lineament density-based risk factor; and [FLULC2020/2030/2040_Risk] indicates future land use and land cover 2020, 2030, and 2040 risk factor value.

FPLR may be subject to change when the new time-oriented data of lineaments, LD maps, and more recent earthquake epicenters are incorporated in future. Therefore, the FPLR maps of three different times such as 2025, 2030, and 2040 in this study, have some limitations only when they are considered in the above-mentioned data layers in future. However, the present research studied the DCEN area to give an idea of how LR will change in future when the future LULC changes in this area. This idea of developing FPLR maps is considered as a pioneering approach for the monitoring of the city sustainability of the DCEN area.

In this stage, for better visualization purposes, a further classification procedure of the natural breaks (Jenks) method with a five-class range (graduated color ramp applied with the blue to red intensification) was implied for each FPLR map. Finally, the prepared FPLR maps of 2025, 2030, and 2040 were masked based on the study area boundary. The entire mapping phases were performed in the multiple software domain, such as Clark Lab’s LCM of TerrSet 18.21, ArcGIS 10.6 including the previously adopted ENVI 5.3, PCI Geomatica 2018, and RockWorks 16 version software, considered for earlier map preparation.

3. Results

3.1. Validation of Future LULCC Prediction

The actual land use of 2018 was compared with the simulated 2018 land use based on the CA-Markov model.

Table 3. Validation of the land use land cover change (LULCC) prediction based on the actual and simulated 2018 LULC image.

Chi-Square Test
LULC Classes	Simulated Land Use in 2018 (O)	Actual Land Use in 2018 (E)	(O − E)2/E
BU	14.94	18.49	0.68
F	38.63	43.73	0.60
AG	15.56	9.61	3.69
RA	0.89	5.17	3.55
WB	5.75	5.42	0.02
BL	24.23	17.58	2.52
Total	100.00	100.00	11.05

Note: χ² = Σ(O − E)²/E = 11.05; df = 5 and χ² _0.05 (5) = 11.07; Data shown in the 2nd and 3rd columns as in O and E represent the percentage share of individual LULC classes used for model validation.

and Figure 4 show the chi-square (χ²) test result for the validation of the model. We hypothesized that the area statistics of the actual and the simulated image were the same.

However, this does not necessarily validate the agreement on the spatial distribution of the LULC classes of the study site. To solve this problem, we performed a more sophisticated Kappa index of agreement between the two images. Moreover, the Kappa coefficient value was measured using the following set of conditions: < 0 = less than chance agreement, 0.01–0.40 = poor agreement, 0.41–0.60 = moderate agreement, 0.61–0.80 = substantial agreement, and 0.81–1.00 = almost perfect agreement. According to Mukherjee et al. 2009 [99], these statistics measure the goodness of fit between the model predictions and reality, which is corrected for accuracy by chance.

The results of the model evaluation performed by the Clark Lab’s LCM module of TerrSet 18.21 version software represents a validation analysis of the agreement/disagreement components (

Table 4. Results of validation analysis of two images (agreement/disagreement components values).

Agreement/Disagreement	Value	Value (%)
Agreement Chance	0.1429	14.29
Agreement Quantity	0.0961	9.61
Agreement Strata	0.0000	0.00
Agreement Gridcell	0.5127	51.27
Disagreement Gridcell	0.1665	16.65
Disagreement Strata	0.0000	0.00
Disagreement Quantity	0.0818	8.18

), which is further partitioned into 0.0818 (error due to quantity/DisagreeQuantity) and 0.1665 (error due to allocation/DisagreeGridcell). Therefore, the data table inferred that the main disagreement between the two maps was due to an allocation error rather than quantity errors between the simulated and actual 2018 LULC images. However, the tabulated chi-square (χ²) value was found to be greater than the calculated one; therefore, we failed to reject the null hypothesis (Table 3). Thus, the CA-Markov model was a good fit to run the future prediction of the LULC for the study site.

Moreover, the overall accuracy of the prediction based on the CA-Markov model could be obtained from the Kno index, which is the standard Kappa index of agreement. The Klocation index validates the simulation to predict the location. All these index of agreement results are shown in

Table 5. Accuracy assessment of the simulated land use image of 2018.

Index	Value
Kno	0.71
Klocation	0.76
KlocationStrata	0.76
Kstandard	0.67

Source: Value calculated by the LCM module of TerrSet Software.

, and the average value is found to be 0.73, which means that the LULC categories of the actual and simulated image were more than 70% similar.

3.2. FLULC Change, Dynamic Degree and Gain and Loss Estimation based on 2018–2040 data

Based on the developed successful simulation of LULC changes in 2018 (Figure 5a), the FLULCC for the different time periods 2025, 2030, and 2040 was generated (Figure 5b–d). The descriptive statistical results (

Table 6. Dynamic degree (DD) (%) estimation between different times based on LULC total area coverage (km²) from 2018 to 2040.

	Total Projected Area Coverage (km2)				DD (%) between Different Times
LULC Classes	2018	2025	2030	2040	2018–2025	2025–2030	2030–2040	2018–2040
BU	58.83	51.06	49.56	48.99	−1.11	−0.59	−0.11	−0.76
F	139.20	122.97	122.98	122.99	−1.67	0.01	0.01	−0.53
AG	30.58	42.30	45.00	46.08	5.47	1.28	0.24	2.30
RA	16.47	1.84	2.28	2.40	−12.69	4.78	0.53	−3.88
WB	55.96	81.28	19.76	79.16	6.46	−0.37	−0.07	1.88
BL	17.26	18.85	18.72	18.68	1.32	−0.14	−0.02	0.37
Total	318.30	318.30	318.30	318.30

Note: BU—built-up area; F—forest area; AG—agricultural area; RA—reconstructed area; WB—water bodies; and BL—barren land. Source: Image statistical results retrieve from ENVI 5.3, LCM of IDRISI-18.21, and ArcGIS 10.6 software based on future projected images of 2025, 2030, and 2040 and DD (%) values calculated by the authors considering Equation (1).

) revealed that from 2018 to 2025, DCEN will have a higher DD in AG (5.47%) than the three other time differences, suggesting faster changes of LULC in the study area (R² = 0.982) (Figure 6c). The BU area shows a negative DD with a −1.11% yr⁻¹ decrease in the time period of 2018–2025 compared to the next three time periods, indicating that the areas of BU will be at risk for future earthquakes.

However, the change rate decreases from −0.59 % yr⁻¹ for the period 2025–2030 to −0.11% yr⁻¹ for the period 2030–2040, respectively, whereas for 2018–2040, the overall change in BU shows a negative trend of development of −0.76% yr⁻¹ (R² = 0.981) (Figure 6a). On the other hand, F also shows a negative change of −1.67% yr⁻¹ in the first phase for the period 2018–2025.

However, Nath et al. [16] reported that after the earthquake, negative changes of DD with −2.03% yr⁻¹ decreased from 2008 to 2018 due to landslides induced by the 8.0 M_w earthquake and two recorded aftershocks, including expansion of the city in the north which resulted in forest loss within the DCEN. The scenario of F shows a positive DD 0.01% yr⁻¹ for 2025–2030, which will remain stable in 2030–2040 and overall will show as negative DD trend with a decrease of −0.53% yr⁻¹ (R² = 0.933) over the period 2018–2040 (Table 6 and Figure 6b).

The FLULC change will occurr in RA with a high negative DD change of −12.69% yr⁻¹ in 2018–2025, and in the next two time periods of 2025–2030 and 2030–2040, positive changes in DD will be observed—4.78% and 0.53% yr⁻¹, respectively. Overall, in the 2018–2040 time period, this trend will negatively impact RA and DD, leading to a change of −3.88% yr⁻¹ (R² = 0.922) in DCEN, which is clearly evident in different time intervals. From 2007 to 2018, it showed a −4.67% yr⁻¹ decrease, whereas in the 2018–2025 and 2025–2030 periods, it decreased −12.69% yr⁻¹ and increased 4.78% yr⁻¹, respectively. By contrast, BL is projected to increase 1.32% yr⁻¹ in 2018–2025 (Table 6). Moreover, in future, WB will change significantly, with an increasing trend of 6.46% yr⁻¹ in 2018–2025 and a decreasing trend of −0.37% and −0.07% yr⁻¹ in 2025–2030 and 2030–2040, respectively. Overall, WB will by 1.88% yr⁻¹ (R² = 0.912) in 2018–2040 (Figure 6e).

With the above changes, BL also shows the same trend of WB, which is an overall change of 0.37% yr⁻¹ (R² = 0.903) in 2018–2040 (Figure 6f). The overall tendencies of WB suggest that if any high-intensity earthquake occurs in this area, WB will probably migrate or the area may be flooded, and also, the frequency of landslides events may increase due to tectonic influence. To crosscheck FLULCC, we performed a change difference of FLULC based on a similar time interval where in each stage, two earlier images and a recent one (for example, 2018 as the earlier ones and 2025 as the recent one) were considered to study changes. This was performed by the minus operator of the Spatial Analyst toolbar of ArcGIS 10.6 software. The change difference maps in different time periods are shown in Figure 7a–d, where a light green color represents ‘no change’ and brick red and blue colors represent ‘negative’ and ‘positive’, respectively.

Further, the future temporal LULC area coverage percentage (%) with the gain/loss (%) estimation was compared for the periods 2018–2040 (

Table 7. Gain/loss (%) estimation between different times based on LULC total area coverage (%) from 2018 to 2040.

LULC Classes	Total Projected Area Coverage (%)					Projected Gain/Loss (%) between Different Times
	2018	2025	2030	2040		2018–2025	2025–2030	2030–2040	2018–2040
BU	18.49	16.04	15.57	15.39		−2.45	−0.47	−0.18	−3.10
F	43.73	38.63	38.64	38.64		−5.10	0.01	0.00	−5.09
AG	9.61	13.29	14.14	14.48		3.68	0.85	0.34	4.87
RA	5.17	0.58	0.72	0.75		−4.59	0.14	0.03	−4.42
WB	17.58	25.54	25.06	24.87		7.96	−0.48	−0.19	7.29
BL	5.42	5.92	5.88	5.87		0.50	−0.04	−0.01	0.45
Total	100.00	100.00	100.00	100.00

Note: BU—built-up area; F—forest area; AG—agricultural area; RA—reconstructed area; WB—water bodies; and BL—barren land. The per year change is not displayed in this table; the calculated result is referred to in the relevant text. The positive value refers to projected gain (%) and the negative value refers to projected loss (%) during different time periods.

) and is shown in Figure 8a, followed by the per year rate of change (%) with trend lines for different periods of FLULC (Figure 8b). To display the trendlines of the future different time periods, a 2nd-order polynomial regression curve (shown in four different colors) was used with regression equations and R² values. This predicted model data help us to achieve a minimum error. Nath et al. [16] earlier reported that in 2007, prior to the earthquake (BEQ) in 2008, the major percentages of LULC in the DCEN were F (39.49%), BL (21.81%), AG (19.80%), and BU (12.74%), whereas in 2008 (AEQ), F (54.90%) increased due to massive plantation activities by the Chinese and local county government.

However, the modeling results show different scenarios with regard to projected changes in LULC in future. The projected BU result shows the coverage area will be 15.39% in 2040 and 16.04% in 2025, both of which are lower than 2018, 18.49%. Overall, the BU data suggest a decline of −3.10% in 2018−2040 (Table 7 and Figure 8a) to −0.14% yr⁻¹, compared to −0.47% (−0.09% yr⁻¹) during the 2025–2030 period and −0.18% (−0.02%) during the 2030–2040 period (Table 7 and Figure 8a,b). By contrast, there is an overall projected loss of −5.09% (−0.23% yr⁻¹) and −4.42% (−0.20% yr⁻¹), respectively, for F and RA during the 2018–2040 period. Compared to the decreasing scenarios, AG, WB, and BL suggest an increase of 4.87% (0.22 yr⁻¹), 7.29% (0.33% yr⁻¹), and 0.45% (0.02% yr⁻¹), respectively (Table 7, Figure 8a,b).

Based on the polynomial regression data fitted to a polynomial equation in three simulated time periods and the overall scenario, the regression equation and R² values show each trend line. During the period 2018–2025, the R² value is 0.234, whereas R²= 0.358 and R² = 0.298 for the periods 2025–2030, and 2030–2040, respectively. However, during the overall 2018–2040 periods, the R² value is 0.251. Our model-based results given in Table 7 and Figure 8b suggest that the 2008 earthquake affected the DCEN area that will be undergoing change which will continue in the future predicted time period up to 2040, unless proper LULC planning and effective measures to minimize the multifarious disastrous risk are not adopted at the policy and decision-making level immediately.

3.3. Analysis of FLULC Transition Probability Matrix (TPM) by Percentage

In reference to Table 2, the 2018 LULC simulated image was generated based on the two available past classified LULC images of 2007 and 2018, adopted by Nath et al. [16]. By considering the 2018 simulated result, attempts were made to classify and generate the FLULC of three time periods, 2025, 2030, and 2040 (Figure 5a–d). The bolded values along the TPM diagonal axis indicate the TP of LULC type will remain unchanged from time (t1) to time (t2), whereas the off-diagonal TPM axes show that a given LULC will undergo a change from one category to another. The derived changes were clearly depicted through the 6*6 FLULC class matrix table, where rows represent the previous LULC categories during time (t1), and the column represents the later LULC classes (t2). For example, the row in 2007 represents actual LULC classes, while the column refers to 2018 LULC classes and thus helps in the prediction of 2018 LULC (Table 2). Similarly, the rows represent the year 2018 and the column represents the 2018 (simulated) image year and predictions for the year 2025 (

Table 8. Markov chain TPM of FLULC classes for the period 2018–2040 by percentage.

Given:	Transition Probabilities of Changing to:
FLULC Classes	Year	BU	F	AG	RA	WB	BL
	2018–2025	0.9763	0.0002	0.0004	0.0009	0.0112	0.0111
BU	2025–2030	0.9888	0.0000	0.0000	0.0000	0.0080	0.0032
	2030–2040	0.9900	0.0000	0.0000	0.0000	0.0072	0.0028
	2018–2025	0.0000	0.9943	0.0020	0.0027	0.0003	0.0007
F	2025–2030	0.0000	1.0000	0.0000	0.0000	0.0000	0.0000
	2030–2040	0.0000	1.0000	0.0000	0.0000	0.0000	0.0000
AG	2018–2025 2025–2030	0.0737 0.0495	0.0000 0.0000	0.8456 0.8928	0.0047 0.0007	0.0000 0.0000	0.0760 0.0571
	2030–2040	0.0355	0.0000	0.9230	0.0005	0.0000	0.0410
	2018–2025	0.3489	0.0000	0.0012	0.6339	0.0000	0.0161
RA	2025–2030	0.2549	0.0000	0.0000	0.7334	0.0000	0.0117
	2030–2040	0.1443	0.0000	0.0000	0.8487	0.0000	0.0069
	2018–2025	0.0000	0.0000	0.0000	0.0010	0.9990	0.0001
WB	2025–2030	0.0000	0.0000	0.0000	0.0000	1.0000	0.0000
	2030–2040	0.0000	0.0000	0.0000	0.0000	1.0000	0.0000
	2018–2025	0.0000	0.0002	0.0033	0.0002	0.0000	0.9963
BL	2025–2030	0.0000	0.0000	0.0000	0.0000	0.0000	1.0000
	2030–2040	0.0000	0.0000	0.0000	0.0000	0.0000	1.0000

).

Furthermore, similar to those, two predicted LULCs of 2030 and 2040 were generated. All these future TPMs were successfully produced in the LCM module of TerrSet Software during the post-classification comparison of the future pairs of image-based results—2018–2025, 2025–2030, and 2030–2040. The results directed at “from-to” change information identifying how many changes and where will occur in the above-mentioned time period are shown in an integrated manner (Table 8). The results highlighting the FLULC changing status in the coming time period show multiple categorical dynamic change conditions.

The FLULC data for the period 2018–2025 show that the BU class has a probability of 1.12% and 1.11% to change into WB and BL, respectively (Table 8), which shows the future disastrous activities in the DCEN area. The highest level of FLULC transformation is observed in RA, where TP 34.89% will transform into another class, such as BU, in 2018–2025, with a smaller transformation of 25.49% and 14.43% into BU during the periods 2025–2030 and 2030–2040, respectively. The RA will remain as RA with a TP of 63.39%, 73.34%, and 84.87%, respectively, during the three time periods. The F, WB, and BL classes show the highest stability to remain a F, WB, and BL with a probability of 100% for each category during the 2025–2030 and 2030–2040 time periods. Similarly, the BU class was also found to be stable with a TP of 97.63%, 98.88%, and 99.00% to remain as the BU class during the three different time periods. The AG class has a TP of 84.56%, 89.28%, and 92.30%, respectively, for the three time periods to remain as AG and a much lower TP of 7.60%, 5.70%, and 4.10% to be converted into BL. This change is obvious and always has a tendency to put pressure on AG, which likely to lead to an increase in the urban and rural population in the DCEN area.

3.4. FPLR Area Identification, Mapping and Pattern of Change Analysis

In the final level analysis, this study further utilized MC-AHP approach to evaluate the FPLR zone for the periods 2025, 2030, and 2040. In a recent study, Nath et al. [16] suggested that the study area is undergoing stress and shows various LR zones with a particular focus on risk with LULC. The LR data of 2018 help to discover the FPLR areas with some limitations because of the changing nature of the lineaments, faults, and future earthquake locations, which will have significant impact on the landscape of the DCEN. By contrast, the structure of the DCEN is dominated by numerous faults, as shown in the present study.

For future landscape planning and sustainability of the DCEN area, it is important to prepare FPLR maps (Figure 9a–c) as well as to identify the areas which will be at risk in the future. These FPLR maps are prepared by considering all environmental parameters with risk priority details based on the methods adopted from Nath et al. [16]. The DCEN area was identified as an earthquake-prone area in the past as well as recently by observing the ongoing stress pattern and frequency of earthquakes in this region. This may also pose a serious threat in the coming decades, too. Therefore, the study area needs emergency attention to assess the FPLR areas at the landscape level. To prepare FPLR maps, multiple landscape factors such as geology, slope, aspect, fault lines, lineament, and its 2018 density (Figure S1a–g), as well as the FLULC maps of 2025, 2030, and 2040 were considered followed by the integrated weighted overlay (IWO) method, which was performed in each FPLR map preparation.

Moreover,

Table 9. FPLR identification and pattern of change analysis based on integrated datasets in different zones. (a) FPLR of 2025; (b) FPLR of 2030; (c) FPLR of 2040.

a.
Total Risk Weight Class Range	Risk Category	Risk Area (km2)	Risk Area (%)	Geological Zone under Risk	Landscape Type under Risk (Priority Based)	Important Locations under Risk
13–18	Very Low	6.53	2.05	Zone- III, IV, II	BU, AG, WB, F, BL	4, 5
18–20	Low	16.25	5.11	Zone- IV, III, II	BU, BL, WB	4, 8, 6, 5, 7
20–23	Medium	116.21	36.51	Zone- IV, III, II, I	BU, F, AG, BL	6, 3, 4, 5, 8, 3, 2
23–25 25–31	High Very High	110.21 69.10	34.62 21.71	Zone-I, IV, II, III Zone- II, I, IV	F, BU, AG, BL WB, AG, F, BU, BL	1, 2, 3, 4, 7, 6, 8 2, 3, 4, 8, 7
Total		318.30	100.00
b.
Total Risk Weight Class Range	Risk Category	Risk Area (km2)	Risk Area (%)	Geological Zone under Risk	Landscape type under Risk (priority based)	Important Locations under Risk
13–18	Very Low	7.28	2.29	Zone- III, IV	BU, AG, WB	4, 5
18–20	Low	17.78	5.59	Zone- IV, III, II	AG, WB, BU, BL	5, 6, 7, 8, 1, 3
20–23	Medium	112.02	35.19	Zone- IV, III, II, I	BU, AG, F, BL	4, 5, 6, 7, 8, 3, 1
23–25 25–31	High Very High	111.12 70.10	34.91 22.02	Zone- I, II, IV, III Zone- II, I, IV	F, AG, BL, WB WB, F, AG, BU, BL	1, 2, 3, 5, 7, 8 3, 2, 1, 8, 7, 6
Total		318.30	100.00
c.
Total Risk Weight Class Range	Risk Category	Risk Area (km2)	Risk Area (%)	Geological Zone under Risk	Landscape Type under Risk (Priority Based)	Important Locations under Risk
13–18	Very Low	7.28	2.29	Zone-III, IV, II	BU, AG, BL	4, 5
18–20	Low	75.10	23.59	Zone- IV, III, II, I	BU, F, AG, WB	4, 5, 6, 7, 1
20–23	Medium	111.38	34.99	Zone- IV, I, II, III	F, AG, BU, BL	8, 1, 4, 5, 6, 7
23–25 25–31	High Very High	54.40 70.14	17.09 22.04	Zone-I, II, IV Zone-II, I, IV	F, AG, BL WB, AG, BU, BL, F	2, 3, 7, 8 3, 2, 1, 8, 7
Total		318.30	100.00

Note: BU—built-up area; F—forest area; AG—agricultural area; RA—reconstruction area; WB—water bodies; BL—barren land; and important locations (1–8) are referred to as: 1 = Yingxiuzhen, 2 = Dabeiping, 3 = Zipingpuzhen, 4 = Puyangzhen, 5 = Xujiazhen, 6 = Dujiangyan, 7 = Juyuanzhen, and 8 = Zhongxingzhen, respectively; Source: Authors: table prepared based on the integrated data overlay techniques, total weighted score, and its classification and area calculations using the raster calculator of ArcGIS 10.6 software, and important location numbering given based on the locations under each risk category.

represents the FPLR area statistics and risk levels prepared based on the developed FPLR maps of 2025 (Table 9a), 2030 (Table 9b), and 2040 (Table 9c). The three final FPLR maps were developed (Figure 9a,c) based on the integrated layer weighted score, which is referred to in this study as AHPA, the particular method of which is discussed in Section 2.3.3. Based on the overall weighted data (OWD), 13 is assigned as the ‘lowest’ and 31 is marked as the ‘highest’ score. The quantiles natural break (Jenks) method with a five-class range has been applied in the ArcGIS 10.6 software environment for classifying the score. The class range was finally designated as 13–18—very low risk, 18–20—low risk, 20–23—medium risk, 23–25—high risk, and 25–31—very high risk.

Table 9a shows a maximum area of 116.21 km² under the medium risk category in 2025 (Figure 9a and Table 9a), which covers 36.51% of the total area of DCEN and will be gradually decreasing from 112.02 to 111.38 km² between 2030 and 2040, respectively, compared to 119.54 km² (37.56%) in 2018 as reported by Nath et al. [16]. Based on FPLR data for the three time periods, a 110.21 km² (34.62%) area was marked as high risk for 2025, which shows a small increase, 111.12 km² (34.91%) for 2030, and a slowdown in 2040 and a share of 54.40 km² (17.09%) of the DCEN compared to its share of 149.23 km² (46.88%) in 2018. By contrast, a 69.10 km² (21.71%) area will be marked as a very high risk area in the year 2025, which will slightly increase to 70.10 km² (22.02%) in the year 2030 and remain almost unchanged at 70.14 km² (22.04%) in 2040. However, in the period 2018, the area was under very high risk for 25.37 km² (7.97%) [16].

In comparison to 2018, the FPLR areas of the very high risk category will grow approximately 2.5 times in 2025–2040. The rest of the areas will be very low to low risk categories in future (Table 9a-c) and (Figure 9a–c), representing 6.28–7.28 km² in 2025–2040 and 16.25–17.78 km² in 2025–2030, further projected to grow 4.5 times on average as a low risk category up to 2040. This result suggests that in future, the risk pattern will gradually decrease, which is expected for any city’s sustainability and free from disaster risk. However, the FPLR table suggests that the study area (DCEN) will not be free from any disaster, as the risk pattern will be progressing until what is reported for the year 2040.

4. Discussion

In this section, we would like to address a few issues that have been observed during FLULC change model predictions in the changing landscape of the DCEN. Firstly, the predicted and observed data are not perfectly matched; therefore, a higher-level accuracy is not achieved. Generally, a higher level of accuracy is expected when the model is run through the LCM module of TerrSet 18.21 version software. We also detected FLULC scenarios, and its trends of change in the percentage and rate of change were observed per year based on different predicted time periods from 2018 to 2040. In connection with this study, the two LULC maps [16] for the years 2007 and 2018 were considered to prepare raster maps and other relevant vector layers to model FLULC and its changing pattern. In addition, we also used all the parameters into FPLR assessment and modeled FPLR scenarios of the DCEN area for the three time periods of 2025, 2030, and 2040, respectively. The CA-Markov model was adopted in this study to predict the FLULC of 2018 (simulated), 2025, 2030, and 2040 for the DCEN area.

Furthermore, the FLULC images were used to generate future potential landscape matrix (FPLM) tables based on the different time periods of 2018–2025, 2025–2030, 2030–2040, and overall 2018–2040. The model and observed data provide over 70% accuracy, which represents a substantial agreement (Section 3.1) of the model validation and its utilizations. Based on the FLULC image output, we produced FLULCC matrix tables of different time periods that show the LULCC transformation from one category to another, as well as its unchanged categories.

In developing and developed countries, the LULC changes and predictive approaches using the CA-Markov model were used across the globe in recent times. However, those studies solely represented the particular LULC changes identified in the selective regions. However, those studies were not focused on the FPLR approaches which were given much less attention and documented less in the past. In recent times, the DCEN area showed changes in LULC for the first time prior to and after the earthquake of 2008 as well as its present status, including the integrated approach to assess the LR pattern [16]. Our study area, SW China of Sichuan province, particularly the DCEN area, frequently experiences earthquakes of different intensity that cause landslides, etc. Figure 2 shows the impact of the earthquakes of varying magnitudes from 12 May 2008 to 8 September 2019. After the 2008 earthquake, Chinese economic development focused deeply on city restructuring, redevelopment, and expansion works due to an inflow migration of people in city areas from neighboring rural areas, thus putting pressure on the landscape of this area. The recent study by Nath et al. [16] shows LULCC in the past, including present LR risk areas and patterns of risk that developed in the earthquake-affected areas. Therefore, based on the earlier study and ongoing changes in stress pattern due to tectonic activities, the present study focused on the FLULCC, its trends of change, and FPLR assessment in the earthquake stress regions. This LULCC modeling and FPLR study is a newly adopted approach in the DCEN area and is particularly needed to check whether the DCEN area will achieve the sustainability goals or not in future along with the projected landscape modifications.

This study first developed and analyzed FLULC mapping and identified the different FLULCC in different time periods, 2018–2025, 2025–2030, 2030–2040, and overall 2018–2040 in the DCEN area of SW China. Here, FLULC gain/loss (%) as well as rate of change per year in the corresponding time frame were computed. Our projected results showed that the BU class has a probability of 1.12% and 1.11% to change into WB and BL, respectively (see Table 8) during the period 2018–2025, which is indicative of future disastrous activities in the DCEN area. Moreover, the highest level of FLULC transformation was observed in RA, with a TP of 34.89% to transform into another class, such as BU in the 2018–2025 time period, and a smaller transformation in 2025–2030 (25.49%) and 2030–2040 (14.43%) time period. The FLULC changes are observed in the central part, north–south (N-S), northeast (NE), and east (E) directions of the study area, which will be reflected in the FLULC maps (Figure 5a–d) as well as in change difference maps (Figure 7a–d) for different time periods. According to the predicted LULC change statistics, further developmental work, including BU area expansion, will continue in future. Moreover, BU also have a high potential risk in future because of an ongoing stress pattern and more earthquake occurrences in the future, which will affect direct modifications of the landscape in the DCEN area.

However, in 2018–2040, the BU data suggest decreasing conditions (−3.10%) of −0.14% yr⁻¹, which will further decerase (−0.47%) at a rate of −0.09% yr⁻¹ during the periods 2025–2030 and −0.18% with a −0.02% yr⁻¹ in the period 2030–2040. The results of the earlier studies [16] suggested that the LULC was changed through modifications caused by natural and anthropogenic activities, such as massive development of infrastructure, BU area expansion, forest loss, etc. Moreover, based on the FLULC maps, the study further examined FLULC transformation during different time periods by getting the change difference matrix tables (showing in percentage) generated during the CA-Markov chain analysis stage through the LCM module of TerrSet software. The individual time-period-based change matrix table is presented in an integrated manner (Table 8).

In addition, this study finally prepared and analyzed FPLR maps for the years 2025, 2030, and 2040 (Figure 9a—c) by applying the IWO technique, with risk priority becoming very high close to adjoining fault lines, and epicenter areas, and lower in the areas farthest from fault lines. The basic assumption followed in this study is similar to the concept developed by Nath et al. [16], which was that the “higher” the LD value, the more the rock fractured in the surface area and vice versa. However, these two parameters have a natural tendency to change in future through the course of time, as is clearly evident from recent studies [16]. Therefore, from this study, it can be inferred that the DCEN area will be at risk in the coming decades with particular activities such as earthquakes, landslides induced by earthquakes, regular landslides, liquefactions induced damage, slope failure, etc.

Therefore, this study further identified the risk according to geological zones in three projected years, 2025, 2030, and 2040. The geological zones II, I, and IV were identified as in the very high risk category (Table 9a), where landscape types WB, AG, F, BU, and BL were found under high stress in the DCEN until 2025. By contrast, in 2030, the risk pattern will change, with a small increase from 69.10% to 70.10%, and the landscape of WB, F, AG, BU, and BL will be under very high stress. However, the risk will remain unchanged until 2040, with the only changes identified at the landscape levels in the order of WB, AG, BU, BL, and F. Among all the landscape risk categories in three years, the medium category risk area will be ranked top, with a share of 36.51% in 2025, and it will gradually decrease in the following two periods, 35.19% in 2030 and 34.99% in 2040. By contrast, the very high, high, and medium category risks will probably be developed in future, as the pace of development will continue closely towards the direction of numerous fault lines, which makes the city become vulnerable.

In this study, the FLULCC modeling results and FPLR areas and probable priority-based risk assessment raised a question on the environmental sustainability of the DCEN area in the coming decades. Thus, a priority-based sustainable solution would be further recommended, as well as checking the unsustainable LULC movements in the study areas. Therefore, by observing the FPLR pattern, it is high time to properly inspect and discuss with the locals what to do and what their concerns are regarding the future landscape change and FPLR, which will be incorporated in the future planning of this area to make the DCEN area sustainable. Moreover, the remote sensing and GIS integrated techniques help to build a model and validate it and are found to be ideal for plausible future LULCC prediction which may further help to identify future probable risk areas at spatiotemporal scale.

5. Conclusions

The present study utilized geospatial data of the two past LULC maps of 2007 and 2018, LR map of 2018, and other associated derived map products used in an integrated manner to develop FLULC maps, future changing trends, and generate FPLR maps of the DCEN. As the DCEN area was previously affected by a high-intensity earthquake (8.0) and its landscape was severely affected afterwards [16], this DCEN area is among the highly seismic-prone areas affected in the SW of China. The present tectonic stress pattern as discussed in the present study is not stable and will be greatly affected in future, as discussed by Nath et al. [16] through LR analysis. In this research, we studied how the landscape of the DCEN area will be changing in future based on the allocated different time periods of 2018–2025, 2025–2030, 2030–2040 and overall projected change in 2018–2040 and quantitatively analyzed and evaluated the changing nature of LULC in the future. The presented results based on the integrated FPLR maps show a future risk pattern observed through the FPLRA along with FPLULC in the DCEN area.

From CA-Markov modeling, changes in LULC will occur in the future in the DCEN area which may be very devastating if the area is affected by high-intensity earthquakes. Therefore, all developmental planning should be focused on finding a nature-based solution that ensures there is less risk and works in conjunction with sustainability goals, with a further suggestion to stop and check developmental activities in the fault and close adjoining areas, and if possible resettle urban dwellers elsewhere, far from the fault lines. The detailed analysis in the present study will be of great use for baseline information to assist city planners, local administrations, and earthquake engineers for future planning and development and technical-know-how about upcoming landscape changes and FLR of DCEN.

An earlier study [16] suggested that the city expansion, LULCT, and its trend continues towards risk-prone areas, which was evident from the recent satellite image of 2018 showing the ongoing tectonic stress level. However, the idea of FLULCC coupled with FPLR assessment used in the present study is a new approach in earthquake-prone areas, a new unique idea that helps to assess risk in a meaningful manner. This study helps us to know the FPLR areas, including their pattern of risk and priority-based landscape risk identification.

Finally, we found that the multitemporal Landsat images of 2007 and 2018 are very useful to generate a LULC simulation image of 2018, followed by a FLULC image for the years 2025, 2030, and 2040 of DCEN using the LCM of TerrSet software. To understand the FLULC, its trends, and future gain and loss estimation of DCEN in the future, we further calculated the matrix of the FLULC in different time domains. However, based on the integrated layers, including the FLULC of 2025, 2030, and 2040, the LD of 2018 (which will change with time) and other ancillary data were recoded with risk priority with the weighted sum method, applied to prepare FPLR maps of DCEN in a consecutive manner for the years 2025, 2030, and 2040. As we observed, DCEN was identified as a risk-prone area in an earlier study [16], something which may continue to be true in future. The FLULC plan for the study area should be made in advance and needs to be incorporated at policy level by city planners and local administrators with a special consideration of the FLULC maps, its change, and FPLR maps. This study of FLULC changes with integrated FPLR mapping techniques was found to be ideal for future sustainable development in the DCEN area, which is under the Belt and Road Corridor (BRC) of P.R.C. China.

The findings of this study might serve as an example for any earthquake-prone city around the world, including the BRC, which will help to exclusively monitor future Silk Road disaster activities. The present study provides future projections about LULC change in the DCEN area, which may be checked in advance for further landscape damage. Geospatial (such as Remote sensing and GIS) integrated software techniques were also found to be ideal, helping to know the FLULC changing status, including in the FPLR areas. Moreover, the present study proved that the integrated techniques have the ability to know FPLR areas and the probable pattern of change in the BRC. Thus, to minimize the FPLR factors of DCEN, it is high time to exchange the results of the FLULC status and its probable change direction, including the FPLR pattern with academics, city planners, national policy makers, and locals to make the DCEN area of the BRC safe, sound, and disaster-proof. Detailed analysis shows a useful future prediction of Silk-Road disaster, which may further help towards achieving the UN 2030 sustainable goals.