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

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


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

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). 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.  30 30 Level-1 GeoTIFF 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]. 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.

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). 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.  [94].
ISPRS Int. J. Geo-Inf. 2020, 9, 134 5 of 25 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]). 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 Mw). 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 Mw 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.

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

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.

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

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

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]: where D represents the DD model, which refers to rate of change; A a is the area in the initial year; A b 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 (%). 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.

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

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 and Figure 4 show the chi-square (χ 2 ) test result for the validation of the model. We hypothesized that the area statistics of the actual and the simulated image were the same. Note:   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, 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.  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), 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 (χ 2 ) 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, 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.   (Table 6) 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 2 = 0.982) (Figure 6c). The BU area shows a negative DD with a −1.11% yr −1 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 -1 for the period 2025-2030 to -0.11% yr -1 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 -1 (R 2 = 0.981) (Figure 6a). On the other hand, F also shows a negative change of -1.67% yr -1 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 -1 decreased from 2008 to 2018 due to landslides induced by the 8.0 Mw earthquake and two However, Nath et al. [16] reported that after the earthquake, negative changes of DD with −2.03% yr −1 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 −1 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 −1 (R 2 = 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.  (Table 6). Moreover, in future, WB will change significantly, with an increasing trend of 6. 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 -1 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 -1 (R 2 = 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 -1 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 -1 , respectively. Overall, in the 2018-2040 time period, this trend will negatively impact RA and DD, leading to a change of -3.88% yr -1 (R 2 = 0.922) in DCEN, which is clearly evident in different time intervals. From 2007 to 2018, it showed a -4.67% yr -1 decrease, whereas in the 2018-2025 and 2025-2030 periods, it decreased -12.69% yr -1 and increased 4.78% yr -1 , respectively. By contrast, BL is projected to increase 1.32% yr -1 in 2018-2025 (Table 6). Moreover, in future, WB will change significantly, with an increasing trend of 6.46% yr -1 in 2018-2025 and a decreasing trend of -0.37% and -0.07% yr -1 in 2025-2030 and 2030-2040, respectively. Overall, WB will by 1.  With the above changes, BL also shows the same trend of WB, which is an overall change of 0.37% yr −1 (R 2 = 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.
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) and is shown in Figure 8a  Further, the future temporal LULC area coverage percentage (%) with the gain/loss (%) estimation was compared for the periods 2018-2040 (Table 7) 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 2 values. This predicted model data help us to achieve a minimum error. Nath et al. [16] (Table 7 and Figure 8a,b). By contrast, there is an overall projected loss of −5.09% (−0.23% yr −1 ) and −4.42% (−0.20% yr −1 ), 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 −1 ), 7.29% (0.33% yr −1 ), and 0.45% (0.02% yr −1 ), respectively (Table 7, Figure 8a,b). Table 7. Gain/loss (%) estimation between different times based on LULC total area coverage (%) from 2018 to 2040.   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 -1 , compared to -0.47% (-0.09% yr -1 ) during the 2025-2030 period and -0.18% (-0.02%) during the 2030-2040 period (Table 7 and Figures 8a,b). By contrast, there is an overall projected loss of -5.09% (-0.23% yr -1 ) and -4.42% (-0.20% yr -1 ), 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 -1 ), 7.29% (0.33% yr -1 ), and 0.45% (0.02% yr -1 ), respectively (Table 7, Figure  8a,b).

LULC Classes Total Projected Area Coverage (%) Projected Gain/Loss (%) between Different Times
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 2 values show each trend line. During the period 2018-2025, the R 2 value is 0.234, whereas R 2 = 0.358 and R 2 = 0.298 for the periods 2025-2030, and 2030-2040, respectively. However, during the overall 2018-2040 periods, the R 2 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.  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.

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

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

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 S1 a-g), as well as the FLULC maps of 2025, 2030, and 2040 were considered Moreover, Table 9 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.    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 2 between 2030 and 2040, respectively, compared to 119.54 km 2 (37.56%) in 2018 as reported by Nath et al. [16]. Based on FPLR data for the three time periods, a 110.21 km 2 (34.62%) area was marked as high risk for 2025, which shows a small increase, 111.12 km 2 (34.91%) for 2030, and a slowdown in 2040 and a share of 54.40 km 2 (17.09%) of the DCEN compared to its share of 149.23 km 2 (46.88%) in 2018. By contrast, a 69.10 km 2 (21.71%) area will be marked as a very high risk area in the year 2025, which will slightly increase to 70.10 km 2 (22.02%) in the year 2030 and remain almost unchanged at 70.14 km 2 (22.04%) in 2040. However, in the period 2018, the area was under very high risk for 25.37 km 2 (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 2 in 2025-2040 and 16.25-17.78 km 2 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.

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 −1 , which will further decerase (−0.47%) at a rate of −0.09% yr −1 during the periods 2025-2030 and −0.18% with a −0.02% yr −1 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.

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.