CA-Markov Chain Analysis of Seasonal Land Surface Temperature and Land Use Land Cover Change Using Optical Multi-Temporal Satellite Data of Faisalabad, Pakistan

Cellular Automata models are used for simulating spatial distributions and Markov Chain models are used for simulating temporal changes. The main aim of this study is to investigate the effect of urban growth on Faisalabad. This research is aimed at predicting seasonal Land-Surface-Temperature (LST) as well as Land-Use and Land-cover (LULC) with a CellularAutomata-Markov-Chain (CA-Markov-Chain). Landsat 5, 7 and 8 data were used for mapping seasonal LULC and LST distributions during the months of May and November for the years 1990, 1998, 2004, 2008, 2013 and 2018. A CA-Markov-Chain was developed for simulating long-term landscape changes at 10-year time steps from 2018 to 2048. Furthermore, surface temperature during summers and winters were predicted well by Urban Index (UI), a non-vegetation index, demonstrating the highest correlation of R2 = 0.8962 and R2 = 0.9212 with respect to retrieved summer and winter surface temperature. Through the CA-Markov Chain analysis, we can expect that high density and low-density residential areas will grow from 22.23 to 24.52 km2 and from 108.53 to 122.61 km2 in 2018 and 2048, as inferred from the changes occurred from 1990 to 2018. Considering UI as the predictor of seasonal LST, we predicted that the summer and winter temperature 24–28 ◦C and 14–16 ◦C and regions would decrease in coverage from 10.75 to 3.14% and from 8.81 to 3.47% between 2018 and 2048, while the summer and winter temperature 35–42 ◦C and winter 26–32 ◦C regions will increase in the proportion covered from 12.69 to 24.17% and 6.75–15.15% of city.


Introduction
Globally, land use monitoring projects have been integral to international climate and environmental science after the start of the Land Use and Land Cover (LULC) project. Urbanization, marked by replacing natural environments with warm absorbing surfaces and houses, results in high urban temperatures relative to rural and sub urban areas [1]. High temperatures occur in central business districts (CBD) and high density residential (HDR) areas [2]. These kinds of temperature changes may have negative environmental and socio-economic effects on built-up areas, including enlarged consumption of heating and air conditioning thus raising energy prices and pollution-related health threats [3,4]. Nonetheless, due to the substitution of impermeable surfaces and buildings as cities expand their coverage and protection is diminished. The association between future temperature and LULC changes predictions needs to be understood for renewable urbanization and making plans.
Surface temperatures and the effect of localized LULC need to be evaluated in order to enhance precision adaptation, prevention, rule design and execution in urban areas.
Several experiments targeted the interaction between LULC trends and LST utilizing remote sensed data without rendering possible predictions [5]. Owing to a mixture of low warm emissivity, weak latent temperature transmission and high heat absorption capacity, this research explained that susceptible surfaces inside urban extents are distinguished by higher temperatures. Conversely, low temperatures are also characteristic of natural landscapes such as wetlands and vegetated areas [6,7]. Research has also looked at natural and long-term historical temperature shifts related to population development [5]. Researchers have used vegetation and urban indices to demonstrate the measurable association between LST and LULC rather than on potential influences from weather trends [8]. For instance [9], describes temperature reduction with the NDVI (Normalized-difference-vegetation-Index) combined with the NDBaI (Normalized-difference-bareness-index) whereas the NDWI (Normalized-difference-wetness-index) is showing rise in temperature with the NDBI (Normalized Difference Built-up Index). The association among LST and variability of indices of LULC is strong, therefore variations in LULC indices such as FVG (Vegetation Fraction) and NDBI have the ability to predict temperature accurately. There is nevertheless a shortage of literature on the usage of LULC indices to predict the potential spatial spread of urban LULC and LST trends.
Despite their success in predicting trends of population development, only one analysis used indices of land cover to estimate projected distribution of soil surface temperature [10]. While [10] the NDVI has been used to estimate residual city normal ecosystems and potential LST values, Normalized difference Vegetation Index is considered to soak at large vegetation fractions, thereby giving a small temperature variation. Previous researches clearly illustrate that Normalized difference Vegetation is a weaker LST predictor than other vegetation and Non-vegetation indices that is, NDBI and ISA (Impervious Surface Areas). In addition [10], single satellite images were used to calculate NDVI to denote the complete season; a technique that is subject to chance, given that a season may differ considerably with a land cover. Therefore, the method has to be changed such as the including seasonal estimates of indices of land cover. In another analysis, calculated LST using a linear regression method on a variety of indices resulting in the Enhanced Build-up and Bareness Index (EBBI), NDBaI, NDBI, NDWI, NDVI, Built-up-Index(BI), Soil Adjusted Vegetation Index (SAVI) and Urban Index (UI) [11]. Thus [10], states that if many variables are included in a linear regression model, the precision of the obtained dependent variable may be affected because of the clatter affected by the collinearity of explanatory factors. Environment predictions are as valuable as they are reliable but hints that a means to estimate LST correctly without errors due to collinearity need to be found first.
To predict changes in LULC and urban expansion, Markov Chain Models have been used [6,12]. For example [9], a Markov Chain analysis for Doha, Qatar, indicated that a 20% rise in built-up areas by 2020 could be expected. Global and local model are widely used to predict temperature excluding urban trend and considering their impact [13,14]. These models are not very fit for understanding regional phenomena as they are at a coarse resolution and therefore need more downscaling [15]. In addition, temperature changes induced by greenhouse gasses are emphasized by global and regional models notwithstanding temperature variations due to the impact of LULC changes. Markov Chain dependent modeling offers an ability to forecast the transition of the environment, offering insights into potential thermal surface characteristics due to vegetation changes [10]. The research is ideal for forecasting variations in temperature at the same spatio-temporal resolution with seasonal variations in land-use land-cover patterns, hence is able to model regional processes such as urban surface dynamics. The Markov Chain model provides tremendous potential for forecasting future LST that needs to be further explored because of earlier achievements in measuring land use land cover shifts relevant effects, usability, flexibility and parsimony. This study is significant for giving direction on how thermal city environments can be projected into future based on the historical patterns of urban development. Much research from various parts of the globe suggest that urbanization contributes to Remote Sens. 2020, 12, 3402 3 of 22 changes in LST but Pakistan still lacks a literature on this issue. The country's meteorological research have largely used large-scale climate models and in situ meteorological data, focused on precipitation and typically targets agricultural impacts [16][17][18]. Remote sensing-based climate research, however, has remained scarce in the region, particularly at the scale of urban microclimates. The same time, urban development estimates based on remote sensing have concentrated mainly on quantifying recent shifts in LULC over the long term [12].
Many of the communities are growing horizontally and vertically in core markets across the world [19]. Pakistan and specially Faisalabad, do not have a suitable LULC management structure. Recently related urbanization to historical first patterns by means of multi-temporal and spectral Landsat datasets then did not forecast potential developments and their effects on the Faisalabad microclimate. Therefore, to the best of our understanding, predictions of potential population development rates, urban land use designs and its impact on LST based on RS (Remote Sensing) have not yet been produced in Pakistan. Therefore, it is important to predict urban development and consequences for the thermal climate of Pakistani cities with information utilizing remote sensing datasets at medium resolution. This has potential to promote local-level adaptation processes, strengthen the decision related to temperature and boost balanced city development that combines possible micro-climate effects of LULC conversions.
Environmental growth is a troubling problem in developing countries like Pakistan. Pakistan has not made a smooth urban transition due to many factors and one of them is the fact that sustainability has not been a consideration. Thus, research on urbanization, urban expansion and the expansion rate is very important for the economic development of Faisalabad and of great interest to the numerous government authorities, town planners and urban scholars around the globe. Regional policy makers and community planners can profit optimally from the work on sustainable expansion. Such research is particularly necessary in order to decide the borders of the new developments and understand urban sprawl across the Faisalabad region. Rapid urbanization has affected the agricultural land and land use patterns generating problems worldwide. Faisalabad Pakistan is not an exception.
The present research identifies seasonal landcover-indices using Landsat 5, 7 and 8 data representing correlations between seasonal (summer and winter) Land Surface Temperature (LST) and seasonal (summer and winter) Land Use Land Cover (LULC) changes in Faisalabad, Punjab, Pakistan from 1990 to 2018. The study also identified the specific indices most suited for prediction of the seasonal distribution of Land Use Land Cover and LST using the CA-Markov-Chain. A further objective was to use seasonal images of summer (May) and winter (November) in land cover indices instead of single season as used in earlier studies when modeling seasonal patterns of land cover as input in the Cellular Automata Markov model.

Study Area
This research was conducted in the industrial city of Faisalabad, Pakistan. It lies from 31 • 25 15" N to 73 • 5 21" E ( Figure 1). The plain fields of Faisalabad are mostly situated on the upper east side of the Punjab, with a height of 604 feet (184 m) above sea level. The city appropriately spreads a territory of roughly 52,142 acres of land while the district approximately covers 1,443,703 acres. In 1901 the population of Faisalabad was 9171 persons, that growing to 70,000 in 1941 and 179,000 in 1951 after partition due to the migration of Muslim refugees. After ten years, the figure had reached 425,000. According to the 1998 census the population had increased to about 2.009 million, exceeding 3.56 million in 2018. The lower Chenab trench is the fundamental water source infrastructure responsible for irrigation across 80% of the cultivated fields. Faisalabad rests on alluvial loess soils with calcareous characteristics, rendering the region extremely productive. The Chenab river flows in the north-west for around 30 km while the Ravi river flows in the south-east at around 40 km from the area. The temperature highs in the city range from 45 • C in summer to 15 • C in winter. The daily mean maximum

Data Acquisition and Image Processing
Landsat 5 Thematic Mapper (TM), 7 Enhanced Thematic Mapper Plus (ETM+) and 8 Operational Land Imager/Thermal Infrared Sensor (OLI/TIRS) medium resolution images were used to calculate LULC and LST. Images of Faisalabad with path/rows of 149/38 were obtained from United States Geological Survey-Center for Earth Resources Observation and Science (USGS-EROS) (http://earthexplorer.usgs.gov/) for two seasons, summer (May) and winter (November) in1990, 1998, 2004, 2008, 2013 and 2018. All these images had less than 10% cloud cover (Table 1). Landsat data were selected because they are easy to access and widely used for classification in LULC and LST analysis. We used cloud free images for analysis as illustrated in Table 1 and Table 2. For atmospheric correction we used the Fast Line-of -Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module in ENVI v_5.4 (Harris Geospatial Solutions, Inc., USA) [22]. Image pre-processing began with layer stacking to create a multispectral image after combining the necessary bands. The images were also calibrated for noise removal [23]. This technique was important to help provide any useful details about the data at the time of data collection, that is, sensor size, sun elevation and atmospheric state.  (Table 1) were used to construct real seasonal predictions. To check the performance of this model, the 10-years interval of satellite images (1998, 2008 and 2018) was changed. Table 1 illustrate that beginning of the summer (all of images from May) and winter (all images from November)

Data Acquisition and Image Processing
Landsat 5 Thematic Mapper (TM), 7 Enhanced Thematic Mapper Plus (ETM+) and 8 Operational Land Imager/Thermal Infrared Sensor (OLI/TIRS) medium resolution images were used to calculate LULC and LST. Images of Faisalabad with path/rows of 149/38 were obtained from United States Geological Survey-Center for Earth Resources Observation and Science (USGS-EROS) (http://earthexplorer.usgs.gov/) for two seasons, summer (May) and winter (November) in1990, 1998, 2004, 2008, 2013 and 2018. All these images had less than 10% cloud cover (Table 1). Landsat data were selected because they are easy to access and widely used for classification in LULC and LST analysis. We used cloud free images for analysis as illustrated in Tables 1 and 2. For atmospheric correction we used the Fast Line-of -Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module in ENVI v_5.4 (Harris Geospatial Solutions, Inc., USA) [22]. Image pre-processing began with layer stacking to create a multispectral image after combining the necessary bands. The images were also calibrated for noise removal [23]. This technique was important to help provide any useful details about the data at the time of data collection, that is, sensor size, sun elevation and atmospheric state. Scanning lines from 2004 and 2008 were eliminated from Landsat 7 (ETM+) in ENVI v_5.4 using Landsat gapfills.sav extension on each band using a triangular approach. After scanning the lines from each band, both bands were opened in Arc Map 10.6 and the bands exported at the same resolution (spatial and radiometric) as Landsat 7 (ETM+). All the maps were geometrically modified using 45 Ground-Control-Points (GCP), toposheet (1:50,000) and aerial-imagery, all GCP obtained on satellite images at the intersection of invariant features and roads. The city boundary was obtained from Urban Unit Remote Sensing and Geographic Information System (GIS) Lab.

LULC Class Abbreviations Description
Central Business District/Industrial area CBD/I. Area Areas with very high density of buildings and a very high proportion of impervious surface that include central business district and industrial areas.
High density residential HDR High density residential areas and areas under residential development (bare or impervious) with low vegetation fraction. Low density residential LDR Established low and medium density residential areas with high vegetation fraction.

Croplands Crops
Areas where intra-urban agriculture is practiced including research sites could be bare in the summer and winter seasons.
Green-spaces G. Spaces Areas covered by grasslands, shrubs and clusters of trees characterized by high vegetation fraction even during the summer and winter season. Water/Wetland Water Areas covered by water bodies or wetlands.
Seasonal satellite images from 1990, 2004 and 2018 were used to develop the model and evaluate the accuracy for predicting land surface variations whereas data from 1990,1998,2004,2008,2013 and 2018 (Table 1) were used to construct real seasonal predictions. To check the performance of this model, the 10-years interval of satellite images (1998, 2008 and 2018) was changed. Table 1 illustrate that beginning of the summer (all of images from May) and winter (all images from November) season in 1990,1998,2004,2008,2013 and 2018. One satellite image from each season was used annually to allow calculation of surface temperature and vegetation & non-vegetation indices for the prediction of seasonal LULC and seasonal LST distributions. The impact of randomness was eliminated by using seasonal averages.

Land Use Land Cover Classification (LULC) Accuracy Assessment and Mapping
LULC maps for the years 1990,1998,2004,2008,2013 and 2018 were obtained from Landsat 5, 7 and 8 reflective bands of 30 m using monitored Support-Vector-Machine (SVM) algorithms. Six LULC categories that is, Central Business District/Industrial Area (CBD/I. Area), High-Density Residential Area (HDR Area), Low-Density Residential Area (LDR Area), Green-Spaces (G. Spaces), Croplands (Crops) and Water/Wetlands (Water) were obtained from every image as indicated in Table 2 [24].
The Support-Vector-Machine algorithm was selected because it does not interfere with the likelihood function and has minimum training data requirements. The Support-Vector-Machine algorithm has established itself in land use landcover classification relative to other classifiers such as the Parallelepiped, Minimum-Distance, Maximum-Likelihood-Classifier (MLC), Artificial-Neural-Network (ANN) and Mahalanobis-Distance classifiers [25,26]. Field observations support training and accuracy assessment in supervised classification. A field survey conducted between March and June 2018 obtained 25 representative GPS points per class except water class. Based on the points were divided into preparation (80%) and testing (20%). Samples of points using a region of interest instead of points increase the reliability of validation results. Specialist expertise and ancillary Land Use Land cover data from previous studies, aerial photos and toposheets were used to establish ground-truth regions for evaluation of classification accuracy. For determining the precision of LULC classifications, the kappa (K) coefficient and overall-accuracy (OA) were used. Post classification [27] variations for every landcover-class from 1990 to 2018 were utilized to measure urban development trends in Faisalabad. Table 3 lists the vegetation and urban indices that were tested for their potential to predict seasonal LST. Table 3 contains of Urban-Indices (UI) calculated using the digital numbers (DN) of specified bands and indices of vegetation, calculated by reflecting the specified bands [11]. As mentioned multiple indices were evaluated for predicting LST. Table 3. Indices of predicting Land-Surface-Temperature (LST) for 2028, 2038 and 2048.

LST Estimation Using Landsat Data
For each year, land surface temperatures were obtained from the Landsat TM and ETM+ (B6) and Landsat OLI/TIRS (B10) thermal bands developed on the dates shown in Table 1. To limit the seasonal influences, the satellite images were taken in the month of May and November. There are two thermal bands (10 and 11) in the Landsat 8 (OLI/TIRS) imagery but only band 10 was used for LST estimation [38]. Land surface temperature recovery includes the transformation of digital numbers (DN) to radiances (L λ ), the measurement of radiance brightness-temperatures (T B ) and the adjustment of emissivity in order to extract surface temperatures from brightness maps [39].
In this research Landsat TM & ETM+ and OLI/TIRS thermal band data from1990, 1998, 2004, 2008, 2013 and 2018 were used to retrieve seasonal (summer and winter) land surface temperatures. Using the Reflectance Toolbox an extension applied to ArcMap 10.6 (ESRI, USA) Equation(1) was used to transform digital numbers (DN) to radiances (L λ ) [23,40]. The tool extracts metadata files from parameters and applies them to the matching thermal data. Brightness temperature was obtained from thermal radiance by means of Equation (2) which is the Landsat channel basic approximation of Planck's blackbody temperature [8,41].
For LST estimation the following series of equations has been used. For each pixel, digital number (DN) was converted into the radiance (L λ ) as follows: Remote Sens. 2020, 12, 3402 where L maxλ is maximum radiance values, L minλ is the minimum radiance values, QCAL max the maximum quantized calibrated pixel value (corresponding to L maxλ in DN (255)), QCAL min is the minimum quantized calibrated pixel value (corresponding to L minλ in DN (01)) respectively; and the metadata of the Landsat images provide their values. Secondly, the L λ values were converted into brightness temperature T B as in: where K 1 and K 2 are constant value and available from the United States Geological Survey (see Table 4) [42]. From every thermal band, we retrieved from spectral radiance and black body the pixel-based land surface emissivity map (ε), as developed [43] and also applied recently [44]. Ultimately, real LST was obtained using Equation (3) after emissivity correction (ε) was applied to the brightness temperature [42].
The λ sign define wavelength of emitted radiance (λ = 11.5 µm), while p is equivalent to 1.438 × 10 −2mk . Using this method, LST was obtained from all the satellite images mentioned in Table 1. The land surface temperature illustrates long term temperature fluctuations in the season and training models in order to predict LST. The average earth surface temperature for 1990,1998,2004,2008,2013 and 2018 was determined by thermal data for the dates seen in Table 1. This was undertaken to test whether the LST were still rising in reaction to urbanization and to determine if potential development could be expected.

Prediction of Temperature by Selecting Variables
A strong correlation between the surface temperature and predictor variables, with no collinearity among the variables is required in seasonal land surface temperature calculation. Degree of correlation with seasonal LST were assessed by indices appearing in Table 3. The seasonal indices having maximum correlation with seasonal LST were chosen to predict seasonal LST using a linear regression model. The association between such variables was also tested to exclude strongly clustered predictors that may trigger collinearity-related errors. Indices strongly correlated with seasonal LST and weakly with each other were used to develop a multi-variate linear model. To evaluate the model's efficiency, we used it to predict the 2018 observed seasonal LST. Mean Absolute Percentage Error (MAPE) -Equation (4) [29] was used to quantify the accuracy.
where T predicted is the model surface temperature and T observed is the real ith pixel of LST reported from Landsat info. The absolute mean percentage is an accuracy prediction metric that represents error in terms of percentage. The model's precision in temperature prediction was measured using the Nash-Sutcliff performance, Root Mean Square Error (RMS), Mean Bias Error and Agreement-Index. The model was applied after accuracy evaluation to predict the seasonal distribution of LST at 10-years Remote Sens. 2020, 12, 3402 9 of 22 intervals for the period from 2028 to 2048. An interval of ten years was selected because the analysis showed noteworthy variations at parallel stages in the time.

Modeling of LULC Changes for 2028, 2038 and 2048
Numerous predictive models have been used to represent the seasonal LULC [45]. The flowchart in Figure 2 describes the technique for predicting potential seasonal LULC and LST distribution from remote sensing data collection utilizing CA-Markov-chain analysis [5,46]

Use of CA-Markov Analysis to Predict Urban Expansion in Faisalabad
IDRISI Selva 17.0 (Clark Labs, USA) is an optimized, certified image processing and Geographical Information Program that offers approximately 300 modules for interactive spatial information analysis and display [47]. The platform provides environmental control, policy support, risk identification, simulation and surface characterization methods. A Markov Chain, is a tool to evaluate adjustments in land use among cycles by a sequence of values that depend on present state [6,48]. The Markov model allows it possible for the structure to progress from the original state i to j over a time period T [9]. The Markov Chain was selected because of its capabilities in predicting time changes LULC [6,45]. In addition, the Markov Chain can predict complex system variations [45]. For other simulations, then, the Markov Chain outputs are used as inputs to generate maps of potential allocation of land usage. Because of its simplicity and parsimony, the Cellular Automata (CA) was nominated to map the spatial dispersal to predicted urban expansion and effect on LST [6,49]. The 2018, 2028, 2038 and 2048 LULC distribution prediction were made using the coupled Cellular-Automata and Markov-Chain the (CA-Markov chain) models in the IDRISI software Selva v_17.0 [1,[50][51][52]. The changes probabilities maps acquired from the Markov Chain analysis were used as inputs to the Cellular Automata (CA) model that maps the future land use landcover distributions. Hence the synthesis of Markov Chain and Cellular Automata exposed spatio-temporal shifts in land use and land cover. We checked their ability to forecast future LULC trends in a dynamic urban environment and applied a CA Markov Chain to make real predictions. We used seasonal LULC shifts between 1990 and 2004 to predict seasonal LULC distribution for 2018 and seasonal LULC transformations between 1998 and 2008 to predict seasonal LULC distributions for 2018 and contrasted both sets of results. The predicted land use landcover was compared with the shape obtained in 2018 using the Support Vector Machine classifications. The forecast performance for 2018 was measured using the Kappa-Agreement-Index (KIA), which measures the degree of consensus between two maps (1990 and 2004, 1998 and 2008) of the same case [51]. Thus, the KIA was used to evaluate performance of the CA-Markov in predictive LULC changes by comparing the LULC map form supervised SVM classification with the modeled map for 2018. After assessing the CA-Markov-Chain model precision, we used LULC seasonal patterns of 1990, 2018 subsequently for predicting seasonal landscape for 2028, 2038 and 2048.

Prediction of LST Distribution in Faisalabad Using CA Markov Chain Study of Land Cover Indices
The Urban Index (UI) was chosen as the best predictor variable of LST distribution as defined in Section 2.6. To avoid restricting the arbitrariness related with one image, an average UI was determined in each of the years 1990, 2004 and 2018 using images collected in May and November ( Table 1)

Statistical Significance of Analyzing Urban Expansion and LST
We checked the statistical significance of predicted seasonal variations in land use land cover and dispersal of LST between 2018 and 2048. The test was applied to coded land use land cover values and to temperature level values derived from 600 points. For each period the LST groups were coded 1-5, whereas the LULC levels were coded 1-6 based on Markov analysis criteria and performance. We used the Shapiro-Wilk statistic to initially test for regularity [47]. Meanwhile, the variations in land use landcover and LST were checked for relevance using the Mann Whitney statistic [29,53]. We checked the Ha hypothesis about the spatial distributions of seasonal LULC and LST is distinct from the alternative Hb hypothesis: in 2018 and 2048 the pairs of LULC and LST were not same [24]. Table 5 presents variations in seasonal LULC distribution recorded at high accuracy using the SVM Algorithm between 1990 and 2018. As shown, the accumulation of built-up areas in the central region and extended from 1990 to 2018. Table 5 describes summer and winter season in two sections. Urban expansion could be investigated from landcover maps and in summer season, CBD/I. area, HDR and LDR was increased from 1.53 km 2 (0.76%) to 6.35 km 2 (3.14%), 8.77 km 2 (4.35%) to 22.23 km 2 (11.28%) and 57.53 km 2 (28.51%) to 108.34 km 2 (53.70%) ( Table 5). Meanwhile, G. Spaces, Crops and water area was decreased from 18.70 km 2 (9.27%) to 16.40 km 2 (8.13%), 114.31 km 2 (56.65%) to 47.92 km 2 (23.75%) and 0.93 km 2 (0.46%) to 0.52 km 2 (0.26%) respectively ( Table 5).The Table 5 winter part indicate that the CBD/I. area, HDR and LDR area was increased from 1.53 km 2 (0.76%) to 6.35 km 2 (3.14%), 8.77 km 2 (4.35%) to 22.23 km 2 (11.28%) and 57.53 km 2 (28.51%) to 108.34 km 2 (53.70%) and G. Spaces, Crops and water area decreased from 13.25 km 2 (6.57%) to 6.74 km 2 (Table 2). For instance, built-up areas expanded exponentially, while Crops and G. Spaces declined from 1990 through 2018.

Seasonal Temperature Changes Observed Satellite Based from 1990 to 2018
Seasonal LST increased between 1990 and 2018 in response to urban development in Faisalabad. as Figure 3 reveals that in 1990, as compared to later years, during the summer season, the region was dominantly covered by temperatures in the range 24-28 • C and 29-30 • C. The 36-42 • C group was dominant in 2018 while lower categories of surface temperature remained in some parts of the main area of the study region. In the central region with CBD / I. area and HDR area, larger temperature increases were found than in the outer areas with crops and G. Spaces.
Seasonal LST increased between 1990 and 2018 in response to urban development in Faisalabad. as Figure 3 reveals that in 1990, as compared to later years, during the summer season, the region was dominantly covered by temperatures in the range 24-28°C and 29-30°C . The 36-42 °C group was dominant in 2018 while lower categories of surface temperature remained in some parts of the main area of the study region. In the central region with CBD / I. area and HDR area, larger temperature increases were found than in the outer areas with crops and G. Spaces.   A summary of the seasonal changes in LST observed in Faisalabad between 1990 and 2018 is given in Table 6. When the city expanded between the seasons the proportion of summer LST in the categories 29-30 °C decreased by around 37% and 17-19 °C categories decreased by 22% in the winter season due to development of CBD I. area and HDR. During the same period, high LST coverage increased nearly 12% in the summer season (35-42 °C ) and nearly 6% in the winter season (26-32 °C ), A summary of the seasonal changes in LST observed in Faisalabad between 1990 and 2018 is given in Table 6. When the city expanded between the seasons the proportion of summer LST in the categories 29-30 • C decreased by around 37% and 17-19 • C categories decreased by 22% in the winter season due to development of CBD I. area and HDR. During the same period, high LST coverage increased nearly 12% in the summer season (35-42 • C) and nearly 6% in the winter season (26-32 • C), suggesting a strong ground warming bias in Faisalabad.  Table 7 indicates a strong correlation between the summer, vegetation and non-vegetation indices and the land surface temperature. The BI, NDWI, MNDWI (Modified Normalized Difference Water Index), UI and NDBI coefficients are at magnitudes greater than 0.5. The other indices showed a lower temperature correlation; for example, the NDBaI had the poorest surface temperature correlation. The other indices displayed negative temperature association; in the summer season, for instance, the EVI (Enhanced Vegetation Index), FVG, NDVI and SAVI had a negative correlation. The Urban Indices had the strongest strong negative correlation with summer FVG (R 2 = −0.7133) but a strong positive correlation with summer temperature (R 2 = 0.9467).  Table 8 indicates a strong positive correlation among LST and BI, UI and NDBI as suggested by coefficients of correlation exceeding 0.5. The other indices showed lowest temperature correlation; for instance, the MNDWI and NDBaI had the weakest correlation and NDVI, NDWI, SAVI, EVI and FVG had negative surface temperature correlation. Urban Indices had the strongest winter temperature correlation (R 2 = 0.9212) and also had a strong negative correlation to winter FVG (R 2 = −0.2353). As the highest correlation was with LST. The Urban Indices have been discovered to be the excellent predictor of urban LST in comparison to different indices.  Figure 5 illustrates regression model of the urban indices for forecasting seasonal LST. The relationship between surface temperature and urban indices for summer and winter season was too strong that is, summer (R 2 = 0.89) and winter (R 2 = 0.92). Therefore, in both seasons, land surface temperature and urban indices had increased and did not affect their relationship due to saturation that disturb indices like NDVI, as the UI continuous to rise with unbounded temperature.

Extraction of LST from Urban Index
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 23 Figure 5. Prediction of seasonal LST from Urban Index by using Linear model.

Use of Urban Indices in Accuracy Assessment of Temperature Retrievals
The linear regression model was evaluated in May and November 2018 on individual Landsat results. These strongly resembled the observed temperature patterns (see Figure 6). Tsat define is the satellite observed temperature and Tmod define modeled derived temperature. The seasonal temperature obtained from the UI was contrasted with that recovered directly collected from Landsat 8's thermal data depending on 600 sampled points throughout the area of study (see Figure 1).

Use of Urban Indices in Accuracy Assessment of Temperature Retrievals
The linear regression model was evaluated in May and November 2018 on individual Landsat results. These strongly resembled the observed temperature patterns (see Figure 6). Tsat define is the satellite observed temperature and Tmod define modeled derived temperature. The seasonal temperature obtained from the UI was contrasted with that recovered directly collected from Landsat 8's thermal data depending on 600 sampled points throughout the area of study (see Figure 1).

Use of Urban Indices in Accuracy Assessment of Temperature Retrievals
The linear regression model was evaluated in May and November 2018 on individual Landsat results. These strongly resembled the observed temperature patterns (see Figure 6). Tsat define is the satellite observed temperature and Tmod define modeled derived temperature. The seasonal temperature obtained from the UI was contrasted with that recovered directly collected from Landsat 8's thermal data depending on 600 sampled points throughout the area of study (see Figure 1).  Visual examination displayed association between the distribution of land use and land cover assessed using the SVM classifier and the distribution of land use and land cover of 2018 projected using the CA-Markov (Figure 7). The model replicates the spatio-temporal dispersal of seasonal forms of LULC as defined by the support vector machine, driven by in-situ-observations. The average Visual examination displayed association between the distribution of land use and land cover assessed using the SVM classifier and the distribution of land use and land cover of 2018 projected using the CA-Markov (Figure 7). The model replicates the spatio-temporal dispersal of seasonal forms of LULC as defined by the support vector machine, driven by in-situ-observations. The average KIA predicted using CA-Markov between the summer and winter LULC and the dispersal calculated by using SVM classifier were 0.88 and 0.85 (Table 9). In the summer, the agreement between the green spaces and the weakest (KIA = 0.78) among the water classes in the two maps was strongest (KIA = 0.89). In the winter, the agreement between the high-and low-density residential area and the weakest (KIA = 0.79) was among the cropland classes in the two maps was strongest (KIA = 0.86).  According to Table 10, residential areas are projected to rise from their present extent through 2048. CBD/I. areas, for example, are projected to expand from 6.89 to 7.92 km 2 between 2028 and 2048. HDR areas are projected to expand from 22.92 to 24.52km 2 . LDR areas are projected to expand from 113,92 to 122,61 km 2 . As the region expands, crops will likely decline in area from 41 According to Table 10, residential areas are projected to rise from their present extent through 2048. CBD/I. areas, for example, are projected to expand from 6.89 to 7.92 km 2 between 2028 and 2048. HDR areas are projected to expand from 22.92 to 24.52 km 2 . LDR areas are projected to expand from 113.92 to 122.61 km 2 . As the region expands, crops will likely decline in area from 41.60 to 33.51 km 2 while G. Spaces fields will decrease from 15.86 to 12.82 km 2 at the same period. In the winter season, CBD/I. area, HDR and LDR area would be increased from 6.89 to 7.92 km 2 , from 229.2 to 24.52 km 2 and from 2028 to 2048 to 113.92 km 2 to 122.61 km 2 . G. spaces and crops regions probably will reduce in size from 12.80 to 9.04 km 2 and 44.30 to 37.11 km 2 . Thus, according to the model projections and expectations, potential development may primarily be marked by growth of HDR at the decrease of water, g. spaces and crops regions.   (Table 11). The model predicted that the temperature ranges from 24 to 28 • C might decrease coverage 13.52-6.32 km 2 whereas the temperature ranges in the 36-42 • C category are expected to rise from 37.87 to 48.67 km 2 from 2028 to 2048. The high temperature category (>26 • C) likely rise at the expense of lower temperature categories as suggested by the winter Figure 9 winter (a-c). Most of the temperature increase will occur in the CBD/I. area, HDR and LDR areas. Most of the area of low temperature (14-16 • C and 17-19 • C) will become high category areas (23-25 • C and 26-32 • C). Many areas particularly in the center and north western part could shift towards higher temperatures (>26 • C) whereas the extent of lower temperature (14-16 • C and 17-19 • C) might decrease as shown in area, could shift towards high temperature (>35 °C ) whereas the extent of low temperature classes (24-28 °C and 29-30 °C) may reduce (Table 11).  The model predicted that the temperature ranges from 24 to 28 °C might decrease coverage 13.52-6.32 km 2 whereas the temperature ranges in the 36-42 °C category are expected to rise from 37.87 to 48.67 km 2 from 2028 to 2048. The high temperature category (> 26 °C) likely rise at the expense of lower temperature categories as suggested by the winter Figure 9 winter (a-c). Most of the temperature increase will occur in the CBD/I. area, HDR and LDR areas. Most of the area of low temperature (14-16 °C and 17-19 °C ) will become high category areas (23-25 °C and 26-32 °C ). Many areas particularly in the center and north western part could shift towards higher temperatures (>26 °C ) whereas the extent of lower temperature (14-16 °C and 17-19 °C) might decrease as shown in Table 11. The model predicted that temperature range of 14-16 °C and 17-19 °C might decrease in coverage from 13.99 to 6.79 km 2 and from 71.63 to 60.83km 2 while the temperature range in the 23-25 o C and 26-32 o C category is expected to rise from 35.83 to 46.63km 2 and from 18.87 to 29.67 km 2 between 2028 and 2048.

Discussion
In this research and Cellular Automata Markov Chain models were developed for prediction of seasonal land use landcover and seasonal LST in Faisalabad, Pakistan. The predictive ability of a number of indices for landcover were assessed in relation spatio-temporal temperature. Among a number of indices such as NDBI, FVG and NDVI, the UI was considered as the prior index for

Discussion
In this research and Cellular Automata Markov Chain models were developed for prediction of seasonal land use landcover and seasonal LST in Faisalabad, Pakistan. The predictive ability of a number of indices for landcover were assessed in relation spatio-temporal temperature. Among a number of indices such as NDBI, FVG and NDVI, the UI was considered as the prior index for predicting LST distribution. Urban expansion and spatiotemporal temperature were predicted by using linear regression model. When the projected variables are not related to each other, it is appropriate to use multiple linear regression [10,54]. Through a diagnosis analysis, all predictors were correlated hence we found that UI provided the most robust temperature explanation. We used the UI as a reference for urbanization and its projected spread to model possible types of LST spread. The model projected temperature with an absolute average error of 1.85 • C by comparing LST computed from thermal band and linear model using UI. UI's success in projecting urban development based heating can be clarified by recent research that have found it to be closely associated with a number of urban expansion metrics [55]. For example [55], observed UI increased with building density and decreased with NDVI in Tokyo Bay.
While the association between temperature and UI has not been verified in earlier research, the strong predictive strength found in this analysis is attributed to elevated temperatures in areas with HDR and less vegetation. [54] also reported strong urban indices in water-intensive and residential parts of Sri Lanka and Colombo. Research have also revealed that the use of water and household energy rises with the intensity of urban weather, thus the strong association between UI and LST [3]. Given the complexity of urban LULC spread, the SVM is a high accuracy classifier tested both in 1990 and 2018. The study of Reference [25] also shows that the SVM classifier can generate high-precision maps. The high accuracy of the map is related to the use of relevant digitized areas (rather than points) as ground data for classification, so its accuracy exceeds the standard 80% [26].
The derived maps of seasonal LULC showed that residential areas grew while vegetation and water cover declined in the same region between 1990 and 2018, which is compatible with previous studies [56]. The CA-Markov-Chain reliably predicts seasonal land use land cover types as showed by the strong consensus between the projected map of the year 2018 and the map prepared from the supervised rating. Based on variations in seasonal LULC between 1990 and 2018, the CA-Markov Chain predicted that unless other measures such as green spaces, cropland and water are implemented and identical trends exist, coverage of built-up area will tend to grow until 2048 at the cost of landcover. This conclusion is aligned with global projections that human growth would continue to grow at the cost of green space resulting in development of built-up areas [6,12]. The lowest temperature level region is predicted to decrease in this summer and winter study, whilst the region protected by warmer categories such as 36-42 • C and 26-32 • C is anticipated to increase. In addition to seasonal changes in LULC distribution, the increasing trends would see residential areas expand to the detriment of green spaces and water. The predicted temperature increases due to changes in urban development demonstrated by the land use land cover are also consistent with earlier research [4].

Conclusions
This study is aimed at exploring seasonal land cover indices and building the CA-Markov-Chain for predicting seasonal distribution of LST and LULC in Faisalabad. We found that the SVM algorithm applied to urban areas classified using Landsat 5, 7 and 8 imagery achieved an overall accuracy above 80%. The model achieved the function of closely matching the spatio-temporal distribution of seasonal types of LULC as defined by the SVM, driven by in situ observations. The average KIA predicted using CA-Markov between the summer and winter LULC and the distribution calculated using SVM classifier was 0.88 and 0.85. The proportion of summer season observed in Faisalabad between 1990 and 2018, LST in the categories 29-30 • C decreased by around 37% and 17-19 • C categories decreased by 22% in the winter season due to the rapid development of the CBD/industrial and residential area. During the same period, high LST coverage increased nearly 12% in the summer season (35-42 • C) and nearly 6% in the winter season (26-32 • C), suggesting a strong biased trend of land heating in Faisalabad. The model is restricted because it only consider the effects of urban development on temperature changes and does not account for other factors. Although, the effective mitigation policies and revised urban expansion strategies have slightly shifted the patterns of LST. Overall, the findings of the present research demonstrate the value of medium-resolution satellite data, in forecasting potential land surface temperatures in urbanized areas. We conclude that unless control measures are taken, the acceleration of urbanization will increase warming and lead to higher temperatures in the future. The results of this study are essential to warn urban planners in order to understand the consequences of expansion on potential temperature changes and thermal comfort of urban residents. Future studies are, however, required to investigate the viability of these approaches and techniques at global and national spatial scales. The CA model was used for the identification of spatial distribution changes and Markov Chain analysis were used for the prediction of temporal resolution. In future researcher and policy makers will use this model to make new policies to control and manage the urban growth.