Simulating Land Cover Changes and Their Impacts on Land Surface Temperature in Dhaka, Bangladesh

: Despite research that has been conducted elsewhere, little is known, to-date, about land cover dynamics and their impacts on land surface temperature (LST) in fast growing mega cities of developing countries. Landsat satellite images of 1989, 1999, and 2009 of Dhaka Metropolitan (DMP) area were used for analysis. This study first identified patterns of land cover different scenarios (e.g., vegetation loss to accommodate urban growth) on LST and consequently to devise appropriate policy measures. forest, gardens, parks and playgrounds, grassland, vegetated lands, agricultural lands, and crop fields. Bare Soil Fallow land, earth and sand land in-fillings, construction sites, developed land, excavation sites, open space, bare soils, and the remaining land cover types.


Introduction
Urban heat island (UHI) is considered as one of the main causes of urban micro-climate warming [1].It is defined as an environmental phenomenon where air and land surface temperatures (LST) of urban areas are higher than those of its surrounding areas [2].UHI is associated with a number of local problems such as biophysical hazards (e.g., heat stress), air pollution and associated public health problems.As a result, the development of strategies to mitigate UHI is now a key policy challenge in order to reduce urban micro-climate warming and to enhance local livability, public health, and well-being [3].This research focuses on the LST component of the UHI effect, and estimates LST in a simulated environment.This is due to the fact that UHI is related to the spatial distribution of LST [1].As a result, LST and UHI are often used interchangeably in this paper.Being an important contributor to the UHI effect, it is, therefore, critical to obtain LST as a first and key step to the UHI analysis; and then to simulate future LST so that policies can be undertaken to reduce the UHI effect.
Research has identified that multiple factors contribute to the generation of UHI (e.g., changes in land use such as urbanization, loss of vegetation and water body, etc.).Urbanization, being the main driver of land cover changes, is considered as one of the most significant factors in this regard [4,5].Land cover changes (e.g., from vegetation to impervious cover such as pavement, roofs, asphalt), again, are the main causes of changing LST because each land cover type possesses unique qualities in terms of energy radiation and absorption.For example, impervious surfaces have low albedos and absorb much of the incoming solar radiation.However, they re-radiate the absorbed solar energy in the form of thermal infrared heat at night [3].As a result, the differences in LST between land cover types vary depending on time in a day (e.g., day and night).In addition, this cross-sectional relationship between land cover types and LST also enabled researchers to investigate the impact of land cover changes on LST over time [6,7].
Studies on the impacts of land cover changes on LST have been intensified recently due to the availability of remote sensing databases of an area for different time periods.However, most of these studies have focused on the past land cover changes and their impacts on LST [3].Although the findings from these studies are important, they lack the ability to simulate the impact of policy interventions on future LST.Given that various simulation techniques are available to model future land cover changes in an area, as a result, it is equally possible to model future LST of that area.Such simulation of LST would also enable to conduct theoretical numerical simulations to quantify the first-order effects of proposed mitigation strategies.However, research on the simulation of LST is fairly limited to date.
Based on the above discussion, this research has two objectives: first, to identify patterns of land cover changes between 1989 and 1999; and between 1999 and 2009, and also investigate their impacts on LST; using Dhaka, Bangladesh as a case; second, to simulate land cover changes for 2019 and 2029, and estimate their impacts on LST in respective periods.

Literature Review
Around 50% (or 3.5 billion) of world population is now living in urban areas, which is expected to be 60% (4.9 billion) by 2030 [8].This increasing trend in global urbanization influenced many researchers to investigate the potential impacts of man-made activities on urban thermal environment such as the LST and UHI effect [9].Research on LST is, therefore, indispensible in order to inform the development of sustainable urban policy.This section reviews literature on the LST component of UHI effect, its character, and key measures.

UHI and Its Key Characteristics
The concept of UHI was described by Luke Howard in 1833 and, since then, research on this topic gradually intensified [10].In a broader geographical context, UHI is a metropolitan area, which is significantly warmer than its surrounding rural areas.However, UHI can also exist within a city region-which implies that a portion of urban area is significantly warmer than its surrounding urban areas [11].The temperature difference is usually larger at night, in winter, and under weaker wind condition compared to daytime, summer, and stronger wind condition, respectively [12].However, a different scenario might also exist depending on spatio-temporal context of an analysis.Ohashi and Kida (2002) have shown that under clear skies and light wind conditions, cities are typically warmer than surrounding rural environments by up to 10 °C [13].Souch and Grimmond (2006) summarized their findings as: the UHI (1) is primarily a nocturnal phenomenon, (2) can occur throughout the year, (3) is dependent on weather conditions such as wind, and (4) generally consists of higher temperatures in the urban core and commercial locations, with lower temperatures in residential and rural sites [14].

Causes and Consequences of UHI
Unplanned and haphazard urbanization coupled with the poor building design are the biggest causes of heat islands in cities. Research has shown that UHI is primarily caused by the built environment in urban areas, in which natural areas are replaced with non-permeable and high temperature surface of concrete and asphalt [15].These surfaces absorb the sun's heat more then they reflect it, causing surface temperature and overall ambient temperature to rise.In addition, tall building and narrow streets can heat the air trapped between them.The UHI effect is the result of a number of factors including but not limited to urban geometry (i.e., size shape, height, and arrangement of buildings), thermal characteristics of urban surfaces, scarce vegetation, extensive use of air conditioning during hot weather, industrial process and automobile use, and the amount of built-up area that exists within a city [16].The built-up area of a city (e.g., building roofs and walls, road, parking lots, industrial area, and commercial area) has less moisture content, and therefore, evaporates less water into the air causing high surface and air temperature.The opposite is true in the case of vegetation in urban areas, which works as a natural cooling system [17].
In summary, there are four major urban elements which may influence the micro-climate of a city.They are built-up area, vegetation, bare soil, and water bodies [18].UHI can impact local weather and climate, by altering local wind patterns, spurring the development of cloud and fog, and influencing the rate of precipitation [11].
The effects of urban heat island are manifold.First, the higher temperature in a city brought about by UHI has an adverse impact on air quality [19].Compared to rural areas, cities experience higher rates of heat related illness and death.The heat island effect is one factor among several that can raise summertime temperature levels that pose a threat to public health [20].

Measurement of LST
Researchers have identified UHI using a number of different measures depending on the type of components of comparison made (e.g., spatial, temporal) in their studies.For example, when comparisons are made between two time periods, UHI was identified by measuring LST differences between the periods of an identical location.This method is often applied to measure the impact of urban growth on UHI effect (e.g., temperature differences between pre-urban and post-urban conditions) [21].Others have identified UHI through measuring rural-urban temperature differences [21,22].Therefore, the comparison is made between geographical units rather than between time periods.However, most of these studies are based on observed temperature of places.
More recently, researcher started utilizing airborne or satellite remote sensing data to derive UHI based on LST [23,24].This technology also allows monitoring land cover changes over time.These two types of opportunities (e.g., derivation of surface temperature and monitoring land cover changes), therefore, enabled researchers to investigate the links between land cover changes and LST changes of a given site between two time periods, i.e., the monitoring of UHI effect due to land cover changes [25].
In the recent past, researchers used the National Oceanic and Atmospheric Administration (NOAA) data to derive LST and consequently to measure UHI for studies conducted at the regional scale [26].However, in recent years, the Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) thermal infrared (TIR) data have often been utilized for smaller scale studies [27,28].Other studies have used mixed type of images.For example, Liu and Zhang (2011) analyzed Landsat TM and ASTER images in order to derive UHI intensity in Hong Kong city [11].
Different types of land cover indices have also been developed in order to investigate the correlations between land cover changes and LST.Amongst the various indices, studies have found that normalized difference vegetation index (NDVI), normalized difference built-up index (NDBI), normalized difference water index (NDWI), and normalized difference bareness index (NDBaI) correlate strongly with LST [5,11,29].NDVI is used worldwide to express the density of vegetation, monitor drought, monitor and predict agricultural production, assist in predicting hazardous fire zones, and map desert encroachment [28].The NDVI process creates a single-band dataset that mainly represents healthy biomass.This index output values between −1.0 and 1.0, where any negative values are mainly generated from clouds, water, and snow, and values near zero are mainly generated from rock and bare soil.Very low values of NDVI (0.1 and below) correspond to barren areas of agriculture, rock, sand, or snow.Moderate values represent parks, shrub, and grassland (0.2 to 0.3), while high values indicate temperate and tropical rainforests (0.6 to 0.8) [30].NDWI implies the water content within vegetation and water state of vegetation [31,32].NDBI is sensitive to built-up area.It has been recently used as an indicator to represent the extent of built-up areas [5,33].NDBaI is useful to classify different barren lands [34].These indices are used to classify different land cover types by setting the appropriate threshold values [5].The relationship between LST and NDVI has been extensively documented in literature [35][36][37].The urban thermal environment is related to the reduction in vegetation coverage [28].NDBI is used for the extraction and mapping of built-up areas [33].NDWI and NDBaI are being analyzed for delineating water content in vegetation and identifying bareness of soil respectively [5,34].The intensity of LST is directly related to the rate of urbanization, land use patterns, and building density.LST is related to patterns of land use/cover changes, e.g., the composition of built-up area, vegetation, and water bodies [5].
Rinner and Hussain showed that LST is extreme where commercial and industrial land uses are located in Toronto, i.e., correlation exists between NDBI and LST [27].Similarly, Chen et al. have shown that LST has become more prominent in areas that experienced rapid urbanization in Guangdong Province, southern China.This study also analyzed dynamics in the spatial distribution of heat islands.The distribution changed from a mixed pattern in 1990 (where bare land, semi-bare land, and land under development were warmer than other surface types) to UHI in 2000 [5].Lo and Quattrochi found that land use changes resulted in significant UHI effect at urban boundary in Atlanta Metropolitan Area in Georgia [6].
Xiong et al. found that high temperature anomalies were closely associated with built-up land, densely populated zones, and heavily industrialized districts.They analyzed Landsat TM/ETM+ images, NDVI and NDBI indices for UHI analysis of Guangzhou in South China [28].Xiao et al. reported that impervious surface is positively correlated with LST in Beijing, China.They used Landsat TM and QuickBird images to analyze the correlations [38].Weng and Yang argued that low vegetation coverage is one of the main reasons for UHI effect in Guangzhou city, China.They generated land cover and LST maps from Landsat TM images [39].
Based on the above review, the following observations can be made in terms of measuring LST.First, although LST is influenced by a number of factors related to urban land cover, most studies have examined just one factor in establishing the correlation.As a result, the influence of other factors has not been disentangled from the correlations presented in these studies.Second, in all cases, researchers have studied only the past or present dataset in measuring LST.Little attempt has been made so far to analyze changes in LST based on future land cover simulations.This article aims to address these gaps in the literature.

Simulation Studies on Land Cover Changes
An urban area is a complex dynamic system.The growth of a city depends on numerous driving factors like social, economic, demographic, environmental, geographic, cultural, and other phenomena.Modeling such dynamic systems is not an easy task [20].Different tools have been developed over the years to underpin the modeling of urban growth and land use changes.Some popular packages include Geomod [40], SLEUTH [41], Land Use Scanner [42], Environment Explorer [43], SAMBA [44], Land Transformation Model [45], and CLUE [46].These tools, again, utilize a number of methods in order to model land cover change such as Markov Chain [47], Cellular Automata [48], Logistic Regression [49], and Artificial Neural Network (ANN) [50].
Each tool has its own advantages and disadvantages.For example, Geomod is designed to simulate only a one-way transition from one land cover category to another land cover category [40].Again, each method has its own strengths and weaknesses.For example, Markov Chain is better applicable when the trend of land cover changes is known.However, this method lacks spatial dependency and spatial distribution [47].The Cellular Automata method models the state of each cell in an array depending on the previous state of the cells within a neighborhood, according to a set transition rules [48].The purpose of this research is not to discuss the strengths and weaknesses of all the models used in the literature as they have been discussed elsewhere [51].This research used the Multi-layer Perceptron (MLP) Neural Network method in order to model/simulate land cover changes [52].
The MLP Neural Network method makes its own decisions about the parameters to be used and how they should be changed to better model the data.It undertakes the classification of remotely sensed imagery through a Multi-Layer Perceptron neural network classifier using the back propagation algorithm.The MLP also performs a non-parametric regression analysis between input variables and one dependent variable with the output containing one output neuron, i.e., the predicted memberships [53].One of its main advantages is that it is distribution-free, i.e., no underlying model is assumed for the multivariate distribution of the class-specific data in feature space [54].Neural networks are non-linear and can be conceived as a complex mathematical function that converts input data (e.g., remotely sensed imagery) to a desired output (e.g., a land cover classification) [20,50,53,54].

Case Study Area
This study was conducted in the context of Dhaka, Bangladesh.Dhaka, the capital of Bangladesh, is one of the fastest growing mega-cities of the world [20].Therefore, although it is likely that such rapid urbanization in Dhaka has a major impact on land cover changes and consequently on the urban micro-climate, little is known about these trends so far in this context.The study area for this research covers the extent of Dhaka Metropolitan (DMP) area (Figure 1).
Dhaka is one of the largest urban agglomerations and the most densely populated mega city in the world.It is the home of more than 16 million people with a total land area of approximately 304.16 km 2 .The population of this city has increased by approximately 11 million in the past two decades [20].Due to this population explosion mainly due to rural-urban migration and partially due to natural growth, Dhaka is expanding both vertically and horizontally.As discussed earlier, these expansions have been identified to be the major contributors of LST increase [20].Existing evidences to date also justify the above argument because an upward shift of temperature has been noted in the last five decades, with an abrupt fluctuation in the minimum and maximum temperature levels (Figure 2).

Data Collection
Landsat satellite images (1989, 1999, and 2009) were downloaded from the official website of US Geological Survey (USGS) and used in order to reach the research objectives [55] (see Table 1 for details).The study area is located in the Landsat path 137 and row 44.The Universal Transverse Mercator (UTM) projection (within Zone 46 North) and the World Geodetic System (WGS)-1984 datum were applied to the images.The pixel sizes of the images were 30 × 30 m.Note that the three images were captured in different time periods; as a result, the atmospheric conditions might be different.However, the images were processed to level-one terrain-corrected (L1T) product.L1T provides systematic radiometric and geometric accuracy by incorporating ground control points while employing a Digital Elevation Model (DEM) for topographic accuracy.Geodetic accuracy of the product depends on the accuracy of the ground control points and the resolution of the DEM used [55].
In addition to the satellite images, a number of land use maps of the study area were collected in order to verify the accuracy of the classified images.The land use maps represented different time periods (e.g., 1987, 2001 and 2010).The maps were collected from the Survey of Bangladesh, which is a national organization responsible for surveying and preparing map products in the country.Therefore, the maps are generally consistent in terms of scale, units, and methods.However, while the existing maps can be valuable qualitative tools, they can result in the error matrix expressing merely differences between the reference data and the map data classification schemes, rather than map error.The scales of the maps were 1:20,000 including the classes-roads, railways, service, hydrography, boundary, cultivated area, low land, village, residential area, under development area, market/commercial area, educational area, industrial area, park and play ground; airport, graveyard, river, and lake.

Methods
The primary objective of this research is to model/simulate future LST.This was done through simulating future land cover changes, because empirical evidence has shown that a strong relationship exists between land cover indices (e.g., NDVI) and LST, as discussed in Section 2.3.Land cover maps of the DMP area were simulated for the period of 2019 and 2029, based on the patterns of land cover changes from 1989 to 2009.Four land cover indices were derived from the simulated land cover for each period including: NDVI, NDWI, NDBI, and NDBaI.The following sub-sections briefly present each of these steps.

Derivation of Land Cover Maps
The acquired satellite images (1989, 1999, and 2009) were classified into four broad land cover types, as shown in Table 2.A supervised classification method was applied as it relies on a priori knowledge of the study area.All the images were analyzed with respect to their spectral and spatial profiles in order to develop training sites [20].A chosen color composite (RGB = 432) was used for digitizing polygons around each training site for similar land cover.Then a unique identifier was assigned to each known land cover type [52].The training sites developed for this research were based on the reference data and ancillary information collected from various sources.When the digitization of the training sites was completed, the statistical characterizations (i.e., signatures) of each land cover class were developed [20].Then, a maximum likelihood classification method was used with equal priori-probability for all classes.The maximum likelihood classifier calculated for each class the probability of the cell belonging to that class given its attribute values [30].Consequently, a filtering technique was applied to generalize the classified land cover images.This post-processing operation replaced the isolated pixels to the most common neighboring class [52].Finally, the generalized images were reclassified to produce the final version of land cover maps for different years (Figure 3).The classified images were then assessed for accuracy based on a random selection of 200 reference pixels for each time period, which were compared against the collected land use maps [52].The overall accuracies of the classified images (1989, 1999, and 2009) were, respectively, found to be 86.48%,90.69%, and 94.83%, with Kappa coefficients of 0.86, 0.91, and 0.95 (Table 3).Note that the Kappa coefficient is a measure of the proportional (or percentage) improvement by the classifier over a purely random assignment to classes [56,57].On the other hand, the user's accuracy measures the proportion of each land cover class, which is correct whereas producer's accuracy measures the proportion of the land base, which is correctly classified.It is observed that the accuracy of the land cover images is increasing over time.The reason is the availability of more detailed and higher resolution reference maps in recent times.

Retrieval of Land Surface Temperature
Surface radiant temperature is derived from geometrically corrected TM and ETM+ thermal infrared (TIR) channel (band 6).Band 6 records the radiation with spectral range in 10.4-12.5 μm from the surface of the earth [55].The geometrically rectified images are free from distortions related to the sensor (e.g., jitter, view angle effects), satellite (e.g., attitude deviations from nominal), and Earth (e.g., rotation, curvature, relief) [55].
Although the impact of the diurnal heating cycle on the LSTs will be an interesting issue to address, there has been no attempt to include it here because TM/ETM+ images do not provide day and night infrared images at the same day.This is why the variability of LST at overpass time in different years is not analyzed.Moreover, because absolute temperatures are not used for the purpose of computation, atmospheric correction was not carried out at this stage.It means no radiometric normalization has been performed.These can be considered as some potential limitations of this research.
The LST was measured from the individual thermal images and were compared between different time periods.Based on the literature, different retrieval methods of brightness temperature from the TM and ETM+ images were applied as discussed below [5].

Retrieval of LST from the Landsat 5 TM Images
Based on Chen et al. (2002), a two step process was followed to derive brightness temperature from the Landsat 5 TM Images in this research [58].In the first step, the digital numbers (DNs) of band 6 were converted to radiation luminance (R TM6 ) using the following formula: (1) where, V represents the DN of band 6, and (3) In the second step, the radiation luminance was converted to at-satellite brightness temperature in Kelvin, T(K), using the following equation: (4) where, K1 = 1260.56K and K2 = 607.66(mW × cm −2 × sr −1 μm −1 ), which are pre-launch calibration constants under an assumption of unity emissivity; b represents effective spectral range, when the sensor's response is much more than 50%, b = 1.239(μm) [55].

Retrieval of LST from the Landsat 7 ETM+ Images
In this research, a two step process was also used to retrieve brightness temperature from the Landsat 7 ETM+ images based on the literature [5,55].In the first step, the DNs of band 6 were converted to radiance based on the following formula (5) where, information can be obtained from the header file of the images, QCALMIN = 1, QCALMAX = 255, QCAL = DN, and LMAX and LMIN (also given in the header file of the images) are the spectral radiances for band 6 at digital numbers 1 and 255 (i.e., QCALMIN and QCALMAX), respectively [55].
In the second step the effective at-satellite temperature of the viewed Earth-atmosphere system, under the assumption of a uniform emissivity, could be obtained by the following equation: (6) where, T is the effective at-satellite brightness temperature in Kelvin; K1 = 666.09(watts/(meter 2 × ster × μm)) and K2 = 1282.71(Kelvin) are calibration constants; and L λ is the spectral radiance in watts/(meter 2 × ster × μm) [55].

Classification of the Heat Zones
The temperature values derived in 'Kelvin (A)' from the above two processes were converted in to 'Degree Celsius (B)' using the following equation: B=A − 273.15 (7) This research found that the temperature values ranged from 15 °C to 36 °C, which were subsequently categorized into six classes (<18 °C, 18 °C to <21 °C, 21 °C to <24 °C, 24 °C to <27 °C, 27 °C to <30 °C, and ≥30 °C) (see, Figure 4).

Derivation of NDVI, NDWI, NDBI, and NDBaI from the Landsat Imagery
As indicated earlier, four land cover indices were derived in this research including NDVI, NDWI, NDBI, and NDBaI.The following formulas, from the literature, were used to derive their scores [5]: , respectively [55].The NDVI map is illustrated in Figure 5.

Simulating Land Cover Maps for 2019 and 2029
As discussed earlier, there are many land use/cover change modeling techniques available (see Section 2.4).The choice of a right simulation technique/model depends upon individual context (e.g., the processes of land transformation), availability of datasets, objective of the research, and the accuracy of the prediction [20,51].In this research, the combination of Markov Chain [59] and MLP modeling techniques [60] were applied to project the future land cover patterns due to their accuracy in outcomes and wide acceptability [61,62].The MLP neural network model takes into account the rate of changes (e.g., expansion) in the built-up areas at the expense of other land cover types (e.g., in this case water body, vegetation, and bare soil).A great advantage of using this model is that it gives the opportunity to model several or even all the transitions at once.A description on the MLP Markov modeling technique can be found in Ahmed and Ahmed (2012), and as a result, is not described here in detail [52].
As the driving force (change to built-up area) for all these transitions are the same here, therefore it is possible to create a common group of explanatory variables that can adequately model all the transitions [20].First, the observed land cover changes (the transitions from other land cover types to built-up areas) were used as the dependent variables; whereas the spatial variables (distance from all land cover types to built-up areas, water body, vegetation, and bare soil; the empirical likelihood transformation) were used as the independent variables.These two types of variables were used to train the MLP Markov model and to produce the transition potential maps.In this case, the MLP running statistics, based on the training and testing root mean square curves, gives an accuracy rate of 90.8%, which is considered to be good for this type of analysis [52,56].The potential maps depict, for each location, the potential it has for each of the modeled transitions.Second, land cover conversions were predicted using the Multi-Objective Land Allocation (MOLA) algorithm [63].The MOLA looks through all conversions to list the host classes that lose some amount of land and the receiver classes that acquire some amount of land from each host [64].The quantities of conversions were determined by the Markovian conversion probabilities [65].After this, a multi-objective allocation was run to allocate land for all receivers of a host class [63,64].The results of the reallocation of each host class were then overlaid to produce the final predicted maps [20,52].

Simulating LST for 2019 and 2029
A three step method was applied in order to simulate LST for 2019 and 2029, as described in this section.Initially, a correlation (bivariate) analysis was conducted between each land cover index (e.g., NDVI) and the LST for each of the time periods (1989, 1999, and 2009) in order to explore the relationships over the periods.The analysis showed that the bivariate correlations are statistically significant and hold over the years.Second, given that the bivariate correlations were significant, a multiple regression analysis using the derived LST for 2009 as the dependent variable and the four land cover indices (NDVI, NDBI, NDWI, and NDBal) for 2009 as independent variable was warranted in order derive a regression equation to be used for projecting the future indices (2019 and 2029).However, a correlation analysis between the indices (e.g., between NDVI and NDBI) showed that the indices were significantly correlated with each other.As a result, the regression analysis was conducted only based on the NDVI in order to avoid the multicollinearity effect from the regression model.The model estimated a constant and a coefficient showing how the NDVI is related with the LST.Third, given that the explanatory power of the model was found to be representative of previous studies on this topic (i.e., explained more than 95% variance in data), the NDVI was, therefore, used to simulate LST for 2019 and 2029.To do this, the NDVI was simulated first for 2019 and 2029 using the Stochastic Markov model [52].
In general, in Markovian processes the future state of a system in time t 2 can be modeled/predicted based on the immediate preceding state; time t 1 .Therefore, the future state can be predicted not based on the past but rather the present [20,47].Initially the NDVI images were reclassified into suitable classes for Markov chain analysis.It produces a transition matrix, a transition areas matrix and a set of conditional probability images by analyzing two qualitative NDVI images from two different dates (1999 and 2009).Each conditional probability image shows the possibility of transitioning to another class [20,52].The next step is to make one single NDVI map for future prediction aggregating all the Markovian conditional probability images.The final predicted land cover map of 2019 will be based on the past ten year's land cover change pattern on the basis of Markov chain analysis.This prediction is performed by a Stochastic Choice decision model [20].Stochastic Choice creates a stochastic NDVI map by evaluating and aggregating the conditional probabilities in which each class can exist at each pixel location against a rectilinear random distribution of probabilities [54].This is how the NVDI maps of 2019 and 2029 were simulated.Finally the simulated NDVI images were used to produce the LST maps of 2019 and 2029 respectively using the regression equation.

Patterns of Land Cover Changes
Figure 6 shows the distribution of land cover in the DMP area for the three time periods considered in this research.Two trends are clearly distinguishable from the figure: (a) built-up area and bare soil increased gradually over the periods; and (b) water body and vegetation declined gradually.More specifically built-up area increased by 88.78% in the past 20 years, from 1989 to 2009.However, two things are not explainable from Figure 6.First, it is unclear whether the built-up areas were only the receiver of lands from other types or donors as well to other types (e.g., converted from built-up to water body).Second, the patterns of conversions like the extent of conversion from one category to the others).In order to understand these relationships, further analyses were conducted (Figures 7  and 8).    Figure 7 illustrates that there was a slight loss in built-up area category in both time periods (e.g., 1989-1999, and 1999-2009).This means that some parts of the previously existed built-up areas were converted into some other land cover classes.However, a substantial change occurred in the bare soil category in both periods.This land cover category both gained and lost lands with almost an equal amount in each period (Figure 7). Figure 8 illustrates the conversion patterns between the categories.It shows that bare soil was the main contributor to forming the built-up areas followed by vegetation and water body (Figure 8a).In contrast, only water body and vegetation were converted into bare soil.Built-up area was not converted into bare soil type at all (Figure 8b).

Change in Land Surface Temperature
Figure 4 shows the patterns of LST changes in the three time periods (1989, 1999, and 2009).The results show that the LSTs ranged from 15.44 °C to 28.12 °C, 18.39 °C to 34.22 °C and 23.03 °C to 35.73 °C.
Figure 9 shows that no areas in DMP experienced a temperature of ≥ 30 °C in 1989.Around 0.35% areas of the DMP experienced this temperature level in 1999, which was found to increase to around 5% in 2009.In 1989, a larger part of the DMP area (74%) fell within the lower temperature zones (<18 °C to < 21 °C) (Figure 9).But in 1999, a majority of the area (91.40%) was found to fall into the mid-temperature zones (21 °C to < 27 °C).This trend continues, and a larger portion of the DMP area (44%) moved into the higher temperature zones (27 °C to <30 °C) in 2009.Another noticeable fact is that no area remained in the lower heat zones in 2009 (Figure 9).Weng (2001) has highlighted that the best way to understand the impact of land cover changes on LST is to investigate the links between the thermal signatures and land cover types [36].Figure 10 summarises the average LST values by land cover type in the three periods.Clearly, the built-up land exhibited the highest LST, followed by bare soil, water body, and vegetation in all three periods.These findings are similar to that reported by Weng (2001), and Weng and Yang (2004) [36,39].The findings signify that built-up areas increase surface temperature by replacing natural vegetation with non-evaporating, non-transpiring surfaces.
In addition, noticeable from Figure 10 is that the LST increased for all land cover types over the periods (e.g., even for vegetation); again suggesting the urban warming effect-otherwise temperature would have not been increased for vegetation.Ren, et al. also reported similar urban micro-climate warming effect in China.Their study has shown that both rural and urban areas experienced a rapid temperature increase in the last two decades [66].However, in comparison with their rate of temperature increase (e.g., 0.32 °C/10 year), the rate of increase in the DMP was higher as shown in Figure 10.This is probably due to the fact that the DMP area experienced a rapid growth.On the other hand, the differences in the temporal resolution of the images might have an impact on this result.However, the results are indicative enough to support the urban area-warming phenomenon in the DMP area.

Relationship between Temperature and Land Cover Indices
Four land cover indices (e.g., NDVI, NDWI, NDBI, and NDBaI) were derived in order to establish quantifiable relationships between LST and the indices.As documented in the literature, a higher level of LST was found to be associated with a lower NDVI in this research (Figures 4 and 5). Figure 11a presents the relationship between LST and NDVI.In 1989, the NDVI ranged between −0.17 and +0.70, which gradually reduced to between −0.72 and +0.29, in 2009.Therefore, it can be said that the NDVI decreased in DMP over time (Figure 5).This research also found that the NDWI and NDBaI values decreased over the periods as well suggesting that both the water content in vegetation and bareness of soil also declined.Figures 11 and 12 show, the correlations and linear regression equations of the LST and the different indices from 1989 to 2009.The results of multiple correlation and regression analyses indicate that LST presents a positive correlation with NDBI (Figure 12a) and negative correlations with NDVI, NDWI, and NDBaI (Figures 11a,b and 12b).

Simulating the Future Land Cover Dynamics
The main objective of this research was to simulate the future LST in the DMP area.In order to reach the objective, the future land covers of the area were simulated first.As mentioned earlier, the MLP Markov Chain model was used to simulate land cover for 2019 and 2029.The simulated land cover maps are shown in Figure 13.Note that 37% of the DMP area was categorized as built-up in 2009.The simulation results suggest that approximately 49% and 57% of the DMP area will be converted into "built-up area" land cover type in 2019 and 2029, respectively (Figure 13).

Simulating the Future LST Maps
The simulated LST maps of 2019 and 2029 are shown in Figure 14.It is found that approximately 56% and 87% areas of DMP will fall in the highest temperature zone (≥30 °C) in 2019 and 2029, respectively (Figure 14).

Discussion
This research models the impacts of land cover changes on land surface temperature (LST) in Dhaka Metropolitan (DMP) area for the period of 2019 and 2029.Numerous studies have investigated the impact of land cover changes on LST [7,11].On the other hand, a handful literature exist that model urban growth for future [66,67].This study is the first of its kind that attempted to combine both techniques and predicted LST based on simulated land cover changes.In addition to the simulated LST, this research also examined the observed relationship between land cover changes and LST in three time periods (1989,1999,2009) for the same area.The observed relationships were found to be consistent with the well-established findings in the literature that: (a) built-up areas exhibited the highest temperature, whereas the lowest temperature was found to be associated with the vegetated/forested areas within DMP area for all three periods [36]; (b) built-up areas increased significantly (89%) over the periods at the expense of vegetation, bare soil and water bodies [66], and as a result, more areas fell within the higher temperature zones over the periods; (c) a gradual increase in temperature was evident for all land cover types in DMP-suggesting the urban micro-climate warming effect in the area [68]; and (d) regression analysis resulted with a positive association between LST and NDBI; and a negative association between LST and other three land cover indices (e.g., NDVI, NDWI, and NDBaI) considered in this research [35,37].Given that the findings from the observed relationships are consistent with the literature, it is believed that the simulated relationships are likely to be defendable because they are based on well-established methods in the literature-MLP Markov.The simulation results from the model were found to be a real concern for the policy makers of DMP area.The findings presented in this paper suggest that there would not be any spatially distributed sporadic UHI within the study area rather the whole DMP area would become an UHI in Bangladesh.Note that, urbanization is the main driver of land cover changes and consequently LST.However, unless undertaking a radical decentralization policy, it is difficult to stop or reverse the urbanization process in DMP because it contains all the major opportunities (e.g., job) and services (e.g., hospital/education) of the country [20].Growth management policies (e.g., green belt) can be implemented that would contain the growth and consequently help reducing UHI effect [69].In addition, policies must not be limited to horizontal growth management only.As discussed in Section 2.2, vertical growth (e.g., building height) can also cause LST.Additional consideration to implement the new urbanism (e.g., green building) concepts in the planning permission (or development assessment) stage of development would also help reducing the LST in the DMP area [70].
Despite the above findings are indicative, they should be read with caution because the results were derived from satellite images with different time periods.Further research should seek to validate the findings reported in this research using more authentic datasets.In addition, although the UHI and LST were often used interchangeably in this paper, they are not the identical concepts rather LST is considered to be just one contributing factor to UHI.The calculation of UHI requires data from a comparative geography (e.g., rural vs. urban), which was not considered in this research.Further research should seek to extend the study area by incorporating the rural/ less urbanized surroundings of DMP area in order to simulate UHI effect in this context.

Conclusions
This is the first study, to the knowledge of the authors that demonstrated a method of deriving future land surface temperature (LST) of urban areas based on the simulation from observed relationships between land cover and LST changes in Dhaka Metropolitan (DMP) area.In doing so, the following conclusions were reached: (1) the amount of built-up area has doubled between 1989 and 2009, and is expected to increase three-fold and four-fold by 2019 and 2029, respectively; (2) a significant negative association between NDVI and LST where as a significant positive association between NDBI and LST exist which suggest that vegetation reduces the LST and built-up areas increases LST.Therefore, future city planning should focus more on urban greening; (3) more DMP areas are gradually shifting towards the highest temperature zone due to the expansion of built-up areas, and if the current trend continues then almost the entire DMP area will be an UHI in 2029.A compact-town like decentralization of urban areas (e.g., satellite-towns) is, therefore, a possible way forward in order to prevent the formation of large-scale urban heat island effect in the future, (4) the study used one land cover index (e.g., NDVI) and applied a simple regression equation for the derivation of future LST.It is possible that the different land cover indices can be used together as independent factors in a multiple regression equation model in order to derive the LST in a more robust way, possibly through undertaking a factor analysis.Future study should seek to incorporate these and improve upon the model presented here.

Figure 3 .
Figure 3. Land cover maps of DMP area.

Figure 6 .
Figure 6.Percentages of presence of land cover types in DMP area.

Figure 7 .
Figure 7. Gains and losses of the land cover types (Unit: % of Area).

Figure 8 .
Figure 8. Contributions to net change (a) in built-up area and (b) in bare soil of DMP area (Unit: % of Area).

Figure 9 .
Figure 9. Changing pattern of heat zones in DMP area.

3 .
Temperature Variations for Different Land Cover Types

Figure 10 .
Figure 10.The variations of mean LST over different land cover types.

Table 1 .
Details of Landsat satellite images.Source: US Geological Survey, 2012.

Table 2 .
Details of the land cover types.

Table 3 .
Accuracy assessments of the land cover types.