Spatial and Temporal Analysis of Surface Urban Heat Island and Thermal Comfort Using Landsat Satellite Images between 1989 and 2019: A Case Study in Tehran

: Mapping and monitoring the spatio-temporal variations of the Surface Urban Heat Island (SUHI) and thermal comfort of metropolitan areas are vital to obtaining the necessary information about the environmental conditions and promoting sustainable cities. As the most populated city of Iran, Tehran has experienced considerable population growth and Land Cover/Land Use (LULC) changes in the last decades, which resulted in several adverse environmental issues. In this study, 68 Landsat-5 and Landsat-8 images, collected from the Google Earth Engine (GEE), were employed to map and monitor the spatio-temporal variations of LULC, SUHI, and thermal comfort of Tehran between 1989 and 2019. In this regard, planar ﬁtting and Gaussian Surface Model (GSM) approaches were employed to map SUHIs and derive the relevant statistical values. Likewise, the thermal comfort of the city was investigated by the Urban Thermal Field Variance Index (UTFVI). The results indicated that the SUHI intensities have generally increased throughout the city by an average value of about 2.02 ◦ C in the past three decades. The most common reasons for this unfavorable increase were the loss of vegetation cover (i.e., 34.72%) and massive urban expansions (i.e., 53.33%). Additionally, the intra-annual investigations in 2019 revealed that summer and winter, with respectively 8.28 ◦ C and 4.37 ◦ C, had the highest and lowest SUHI magnitudes. Furthermore, the decadal UTFVI maps revealed notable thermal comfort degradation of Tehran, by which in 2019, approximately 52.35% of the city was identiﬁed as the region with the worst environmental condition, of which 59.94% was related to human residents. Additionally, the relationships between various air pollutants and SUHI intensities were appraised, suggesting positive relationships (i.e., ranging between 0.23 and 0.43) that can be used for establishing possible two-way mitigations strategies. This study provided analyses of spatio-temporal monitoring of SUHI and UTFVI throughout Tehran that urban managers and policymakers can consider for adaption and sustainable development.


Introduction
An Urban Heat Island (UHI; see Table A1 for the nomenclature) is a phenomenon in which urban areas experience higher temperatures than surrounding rural areas. This phenomenon is probably one of the most important indicators of environmental conditions in metropolitan areas [1,2]. In particular, UHI adversely affects numerous socio-economic and environmental factors, including urban climate [3], vegetation growth [4], quality of drinking water [5], rain intensity [6], air pollutant concentration [7], human health [8], and energy consumption [9]. The continuous conversions of natural land covers to urban impervious surfaces alter the land surface energy processes and the thermodynamic properties of the surface, increasing the urban temperature and forming greater UHIs [2,10,11].
According to the United Nations report, the urbanization rate has grown from 29% to 55% from 1950 to 2018, and it is predicted that by 2050, 68% of the world's population will indwell in urban areas [12]. This requires rapid urban expansions and developments, which accelerate the temperature rise in metropolitan areas and significantly intensify the UHI formation and intensity [13][14][15]. Therefore, it is of significant importance to monitor and study the UHI patterns in metropolitan areas to adopt effective mitigation strategies and promote sustainable cities [2,16].
Conventionally, in situ observations of air temperature data from ground stations or radiosondes, recorded in both urban and rural areas, were employed to calculate UHI (i.e., at urban canopy layer or urban boundary layer) intensities for further analyses [17,18]. Although in situ observations provide the most precise data, their limited number and spatial discontinuity restricted their usage for efficient UHI mapping [17,19]. Consequently, it is appealing to incorporate remote sensing data, with spatial continuity, broad coverage, and frequent data acquisition, to investigate Surface UHI (SUHI) variability in metropolitan areas [20,21]. Notably, the possibility of recording thermal emissions of the Earth's surface through Thermal Infrared (TIR) remote sensing created an exceptional opportunity to extract Land Surface Temperature (LST) products to study SUHI patterns allowing for a more profound understanding of thermal spatial pattern and the impact of surface characteristics of SUHI [22]. Additionally, the LST products allow calculating different thermal comfort indices, such as Urban Thermal Field Variance Index (UTFVI), a broadly used thermal index for assessing the environmental condition of metropolitan areas, to manifest the impact of SUHI intensities [23][24][25][26].
Up to now, many scholars have dedicated research studies to derive SUHI patterns/maps from remote sensing LST data for further spatial analyses. For instance, Guha et al. [23] employed Landsat-8 data to map the SUHI patterns and UTFVI maps in Naples and Florence, Italy. Their results revealed that more than 75% of SUHIs were formed within Bare Land (BL) and Built-Up (BU) areas, which were also demarcated as ecologically stressed zones. In another study, Bokaie et al. [27] employed Landsat Thematic Mapper (TM) data to map the SUHIs of Tehran in 2010 and investigate its relationship with the Land Use/Land Cover (LULC) map and Normalized Difference Vegetation Index (NDVI) image. They reported full compliance between average LST values and LULC classes and a moderate negative correlation between LST and NDVI values, which was also in accordance with other studies [28]. Likewise, several other scholars incorporated multi-temporal remote sensing data to map the spatio-temporal variability of SUHI patterns [29][30][31]. For example, de Faria Peres et al. [32] explored the trend of SUHI evolution over 30 years and compared the results with LULC maps. The results suggested that the main reason for the 2 • C rise of SUHI intensity in Rio de Janeiro was associated with urban expansion due to the significant growth of LST in urban areas. Additionally, Nadizadeh Shorabeh et al. [33] employed five Landsat images between 1985 and 2017 to study the SUHI variations in Tehran. Later, they applied the Cellular Automata-Markov (CA-M) and Artificial Neural Networks (ANN) model to predict the LULC of 2033 to model the future surface SUHI intensity.
Tehran is the largest and most populated metropolitan in Iran, and as the central hub (i.e., political, economic, social) of the country, it has experienced enormous population growth and extensive urbanization [34]. Several studies were carried out to study and monitor SUHI and LST variations throughout the city [35][36][37][38][39][40][41]. However, the SUHIs were still extracted by a single image in these studies, so that they could not be considered as a thorough description of annual or seasonal SUHI. This is because Utilizing timeseries remote sensing images produces a more detailed and persuasive understanding of the complexity of SUHI in comparison with analyzing this phenomenon with limited images [42,43]. Moreover, the thermal environmental condition of Tehran has not been analyzed in previous studies. To the best of our knowledge, no comprehensive study was dedicated to investigating three decades of SUHI and UTFVI patterns in Tehran through time-series data. Furthermore, Tehran is suffering from severe air pollution [44], and thus, it is required to appraise the relationship between air pollutants and SUHI intensities, which has not been conducted in Tehran. In fact, the contradictory reports of the relationship between air pollutants and SUHI intensities in different locations necessitate performing these analyses for Tehran [45][46][47][48]. These investigations would provide profound information about the environmental condition of Tehran, leading to effective decision-making for a sustainable city.
Considering the foregoing, this paper aims to extend previous studies and provide relevant information from new aspects by investigating the spatio-temporal variability of SUHI and thermal comfort and appraising the relationship of SUHI intensities and air pollutant concentrations in Tehran. Specifically, the present study follows three objectives: (1) Investigating the SUHI changes over the past three decades and examining its intraannual variations, providing the SUHI magnitudes and footprints; (2) exploring the spatial changes of the environmental condition of Tehran over the last three decades using the UTFVI; and (3) identifying the relationship between SUHI intensities and different air pollutants concentration for Tehran.

Study Area and Data
This section covers three subsections, namely the study area, satellite remote sensing data, and ground data.

Study Area
Tehran, covering an area of about 730 km 2 , is the capital of Iran, located between 34 • 35 -35 • 50 N and 51 • 02 -51 • 36 E in the northern part of Iran ( Figure 1). Tehran is surrounded by the Alborz Mountains to its north and the country's central desert to its south, with a significant impact on the formation of a semi-arid climate. The climate is mild in spring and autumn, hot in summer, and cold in winter, especially at night. With a high elevation range between 900 and 1800 m above sea level, Tehran includes diverse annual temperatures, reaching 42 • C in July [49]. Tehran is the central economic, political, and recreational hub of Iran that attracts the most internal immigrants. Being the first immigrant destination, the city population has grown from 6 million in 1985 to over 8.5 million in 2017 [33]. Consequently, the city has undergone significant urbanization and urban expansion, leading to several environmental issues, such as air pollution, heavy traffic, and high energy consumption demands [50,51]. Furthermore, the widespread conversion of natural land covers and Green Space (GS) areas to BU and impervious surface areas, such as altering the thermodynamic characteristics of the surface [10], resulted in an increase in urban temperature and formation of SUHIs with higher intensities [52]. Therefore, profound research should be dedicated to studying and monitoring long-term SUHI patterns and the thermal comfort of Tehran to promote a sustainable city through adaptation and mitigation strategies.

ployed.
In total, 47 Landsat-5 scenes (i.e., Path 164 and Row 35), with less than 25% cloud coverage, between 1989 and 2012 were used for LULC mapping and generating SUHI and UTFVI maps in 1989, 1999. Additionally, 21 Landsat-8 scenes (i.e., Path 164 and Row 35) with less than 25% cloud coverage were also used to derive LST images and map SUHI and UTFVI in 2019 and to produce LULC maps between 2013 and 2019. Table A2  provides brief specifications of the satellite images used in this study, and the ImageCollection IDs are provided in Table A3.

Ground Data
In the present study, field observations from pollutant monitoring stations were employed to investigate the relationships between SUHI intensity and different air pollutants. The hourly measured variables from 20 stations (see Figure 1) were provided by the Tehran Air Quality Control Company. In this regard, the concentration values of air pollutant variables, including NO2, O3, and PM with aerodynamic diameters of 2.5 and 10 microns for the year 2019, were downloaded from (https://air.tehran.ir accessed on 10 January 2020). Downscaled daily average concentration values were utilized for further analyses. Table 1

Satellite Data
TIR and optical (i.e., visible, infrared, and shortwave infrared) bands of Landsat-5 and Landsat-8 satellite datasets were utilized for LULC, SUHI, and UTFVI mapping. Landsat-5 and Landsat-8 are the fifth and eighth satellites of the Landsat Program, which were launched as a joint cooperation between the United States Geological Survey (USGS) and the National Aeronautics and Space Administration (NASA).
Time series data from 1989, 1999, 2009, and 2019 were utilized for a thorough analysis of SUHI and UTFVI throughout the year. These satellite data were collected from the Google Earth Engine (GEE) cloud computing platform, which provides efficient and rapid data access for geospatial data processing and prototyping [53,54]. Here, the surface reflectance of Landsat-5 and Landsat-8 data within GEE through ImageCollection IDs of (LANDSAT/LT05/C01/T1_SR) and (LANDSAT/LC08/C01/T1_SR), respectively, were employed.
In total, 47 Landsat-5 scenes (i.e., Path 164 and Row 35), with less than 25% cloud coverage, between 1989 and 2012 were used for LULC mapping and generating SUHI and UTFVI maps in 1989, 1999, and 2009. Additionally, 21 Landsat-8 scenes (i.e., Path 164 and Row 35) with less than 25% cloud coverage were also used to derive LST images and map SUHI and UTFVI in 2019 and to produce LULC maps between 2013 and 2019. Table A2 provides brief specifications of the satellite images used in this study, and the ImageCollection IDs are provided in Table A3.

Ground Data
In the present study, field observations from pollutant monitoring stations were employed to investigate the relationships between SUHI intensity and different air pollutants. The hourly measured variables from 20 stations (see Figure 1) Figure 2 presents the flowchart of the proposed method for monitoring and mapping the SUHI and UTFVI patterns. The proposed methodology includes four subsections, which are described separately in detail. In this regard, first, the satellite data preprocessing is explained. The LULC and SUHI mapping are described in the two next subsections, followed by the UTFVI explanations in the last subsection.   Figure 2 presents the flowchart of the proposed method for monitoring and mapping the SUHI and UTFVI patterns. The proposed methodology includes four subsections, which are described separately in detail. In this regard, first, the satellite data preprocessing is explained. The LULC and SUHI mapping are described in the two next subsections, followed by the UTFVI explanations in the last subsection.

Satellite Data Preprocessing
As mentioned earlier, Landsat-5 and Landsat-8 surface reflectance datasets were employed. Several preprocessing steps were initially applied to Landsat series satellite images by the GEE developers. Regarding the initial preprocessing, the Landsat Ecosystem Distribution Adaptive Processing System (LEDAPS) and the Land Surface Reflectance Code (LaSRC) algorithms were respectively applied to Landsat-5 and Landsat-8 images

Satellite Data Preprocessing
As mentioned earlier, Landsat-5 and Landsat-8 surface reflectance datasets were employed. Several preprocessing steps were initially applied to Landsat series satellite images by the GEE developers. Regarding the initial preprocessing, the Landsat Ecosystem Distribution Adaptive Processing System (LEDAPS) and the Land Surface Reflectance Code (LaSRC) algorithms were respectively applied to Landsat-5 and Landsat-8 images for atmospheric correction. Afterward, the C Function of Mask (CF Mask) was performed to identify cloud, shadow, water, and snow masks for each pixel [55].
In this study, only cloud-free satellite data over the study area were used for further processing. The surface reflectance optical data were directly used to produce LULC maps and NDVI images of the study area, while other preprocessing steps, including surface emissivity estimation (Equations (1) and (2)) [56] and emissivity correction to derive LST images, were applied. It is worth noting that the emissivity correction increases the reliability of the LST data [57,58] and, thus, improves the reliability of SUHI and UTFVI analyses.
where ε and P v are surface emissivity and vegetation cover density, respectively, and NDVI values were calculated based on Near Infrared (ρ NIR ) and Red (ρ Red ) surface reflectance bands. Subsequently, estimated surface emissivity, along with brightness temperature images, are employed to calculate the LST images by Equation (3).
where L λ is brightness temperature in kelvin, and λ is the effective wavelength of the emitted radiance. Moreover, α is calculated based on h.c/k, in which k is the Stephan Boltzmann constant (1.38 × 10 −23 J/K), h is the Planck's constant (6.626 × 10 −34 J.s), and c is the velocity of light (2.998 × 10 8 m/s).

Land Cover/Land Use (LULC) Mapping
Optical surface reflectance images from both Landsat-5 and Landsat-8 satellites were employed to generate four class LULC maps with 30 m spatial resolution from 1989 to 2019. The produced LULC maps were used to investigate the urbanization trend during the last three decades and also to examine their relationship with the spatio-temporal patterns of SUHI and UTFVI maps for the years 1989, 1999, 2009, and 2019. The Support Vector Machine (SVM) classifier with the Radial Basis Function (RBF) kernel implemented within GEE was utilized to produce LULC maps. SVM is a non-parametric classification algorithm, which is based on fitting optimal separating hyperplanes between different classes by focusing on the training samples that lie at the edge of the distribution of the classes in the features space, called support vectors [59]. It is worth noting that the RBF kernel was used since it provides accurate LULC maps [60], and its necessary tuning parameters were determined based on several trial and error attempts to achieve satisfactory classification results. The generated LULC maps include four classes, namely Water Body (WB, including artificial lakes), Green Space (GS, including urban parks, urban forests, trees, croplands, and grasslands), Bare Land (BL, including unused lands with soil cover), and Built-Up (BU, including human residents, industrial sites, and urban infrastructure). The required reference samples to support the supervised classification tasks were collected by precise visual interpretation of each image, and then they were randomly split into two halves of training and test samples [61]. On average, 248 polygons with an area of about 4.7 km 2 , including four LULC classes, were collected for each period, and Table A4 provides the number of polygons and area of each land cover for the years 1989, 1999, 2009, and 2019.

Surface Urban Heat Island (SUHI) Mapping
Generally, SUHI mapping methods are categorized into three groups, namely (1) using LST as a proxy of SUHI, (2) LST differences between urban and the surroundings [13,62] (i.e., reference areas), and (3) statistical methods [20]. The first approach uses LST values to investigate the SUHI variability since the SUHI manifests itself in hotspot forms or high LST values in comparison to the surrounding environment [63]. The main limitation of these approaches is that they avoid the measurement of SUHI intensities and make SUHI comparisons more challenging [20]. In the second category, the LST differences between urban areas and mean LST values of reference areas (e.g., rural, suburban, water, and vegetation) are employed to map the SUHI [38,64]. Although these methods are simple to implement and provide SUHI intensities, they commonly suffer from uncertainties associated with urban and reference areas delineation [65]. The third category utilized statistical algorithms, such as Gaussian Surface Modeling (GSM) and linear regression functions, to generate SUHI maps [66,67]. These approaches are not affected by the biases caused by reference area selection and, thus, facilitate the SUHI comparisons [20].
In this study, the planar fitting and GSM were implemented to map the SUHI areas of Tehran. The GSM model enables a better understanding of the SUHI patterns as it provides the magnitude, spatial extent, and central location of the SUHI. Furthermore, this method allows quantitative comparisons of SUHI maps over timescales for a single city [68]. The GSM was first developed by [69], and then was employed in other studies to delineate SUHI intensity and its spatial extents. A stable workflow to implement the GSM model for SUHI mapping is thoroughly described in [70]; however, brief step-wise explanations are provided in the following. For this task, the first step is to compute the SUHI map based on Equation (4). To this end, the urban areas (BU pixels from the LULC maps) were temporarily masked from the LST images (T(x, y); x and y are the pixel coordinates), and then a planar surface was fitted to the existing non-urban LST values to determine T 0 , a 1 , and a 2 . This part attempts to model the spatial pattern of LST in the study area with the assumption of no BU areas. The next step is to subtract the fitted surface from the actual LST with all pixels to extract the SUHI patterns (Equation (4)). Afterward, the GSM, as shown in Equation (5), is fitted to the derived SUHI images using the least-square method. These approaches provide the SUHI characteristics, including SUHI magnitude (a 0 ), the central location of the primary heat island (x 0 and y 0 ), its orientation (ϕ), and its spatial extent (a x and a y ). Furthermore, as an automatic workflow, it eliminated the uncertainties associated with reference area selection [20].

Urban Thermal Filed Variance Index (UTFVI) Mapping
UTFVI, which is calculated based on Equation (6) [71], is an extensively used indicator that describes the environmental condition and quality of urban health through the measure of thermal comfort presence in the environment [25]. Although there are different indices for the evaluation of thermal comfort, they are not applicable in some regions due to the lack of data [72][73][74]. UTFVI notably causes adverse impacts on different socio-economic and environmental issues, including humidity, air quality, indirect economic loss, reduced living comfort, and increased mortality rate [75]. Accordingly, the examination of the UTFVI patterns is attracting more interest in the scientific community since it could provide beneficial information of thermal comfort status, ensuring sustainable development [25]. The calculated UTFVI values were then classified into six classes (see Table 2) for better visualization and to provide explicit environmental condition patterns throughout the city [23].  For a more profound comparison between SUHI variations and LULC changes, the SUHI maps were classified into five categories, provided in Table 3. Table 3. Surface Urban Heat Island (SUHI) ranges to define categorial SUHI. Figure 3a-d illustrates that the city has undergone a substantial LULC change, mainly in terms of urban expansion and the conversion of natural land covers into urban impervious surfaces and BU. In particular, the produced LULC maps, with overall accuracies that varied between 89.14% and 93.61% (see Figure A1), indicate that the BU areas covered 329.49 km 2 , 387.24 km 2 , 446.93 km 2 , and 505.17 km 2 of the city in 1989, 1999, 2009, and 2019, respectively. This specifies an approximately 53.33% growth in BU areas over the last three decades. Additionally, the time-series LULC maps generated from 25 Landsat-5 and Landsat-8 images between 1989 and 2019 revealed a gradual upward trend of urbanization with a rate of 5.86 (±2.42) km 2 per year in Tehran. According to Figure 3a-d, considerable urban expansions occurred in the western and northern parts of the city due to the construction and development of residential, industrial, and commercial infrastructures. For instance, in 1989, most of the western parts of Tehran were covered by BL and GS, constantly changing to BU by 2019. Moreover, a notable part of GS areas in 1989, especially in the northern parts, were lost and superseded by BU.
Likewise, Figure 3e-h presents the SUHI maps of Tehran between 1989 and 2019. As is clear, the spatial patterns of SUHIs have changed over the past three decades. Visually, the regions in which the SUHI intensities were increased are mostly demarcated in locations where BU areas were developed or BL areas exist. For example, the northern part of the city experienced SUHI intensification that is mainly associated with the drastic reduction of GS areas (see Figure 3a-d), which had eliminated the cooling effect of this land cover. Similarly, the western regions, which have been the center of urban development for the past decades, encountered a significant increase in SUHI intensity. Considering LULC maps (Figure 3a-d), the SUHI intensifications could be partly identified as the consequences of GS loss across Tehran. In addition to SUHI intensity changes over the regions associated with LULC transitions, the city center, which was almost covered by BU in the study period, also experienced a moderate aggravation of this phenomenon. Despite spatial changes of SUHI intensities throughout Tehran, further statistical analysis revealed that the average SUHI intensity in Tehran increased from 2.05 • C to 4.07 • C from 1989 to 2019. This also confirms the necessity of devoting more research studies to examine the SUHI patterns and intensities for further adaption.
The UTFVI maps were also derived from LST images and then were classified into six classes (see Table 2) for better visualization and to provide an explicit spatiotemporal overview of the environmental condition and thermal comfort throughout Tehran (Figure 3i-l). Generally, major BU areas were distributed over two extreme categories of EC and WTC, while small BU patches were demarcated in GC and NC. As in the SUHI maps and urbanization patterns, the western part of Tehran had a higher rate of UTFVI alterations, in which the WTC has spatially expanded except for small patches, including the Chitgar Lake and the surrounding areas. The center and south of the city have also experienced thermal comfort degradation since the WTC patches have enlarged. UTFVI maps of 1989-2009 illustrate that over 55.45% of Tehran was identified as EC, while in 2019, 52.35% and 37.11% were classified as WTC and EC, respectively. The results manifest a 16.37% and 18.34% increase and decrease in the areas of WTC and EC in 2019, respectively, compared to the same classes in 2009.

Relationship between SUHI and LULC
For a more profound comparison between SUHI variations and LULC changes, the SUHI maps were classified into five categories, provided in Table 3. The classified SUHI maps, along with the corresponding pie charts showing the share of each SUHI class, are presented in Figure 4. It is evident that over the past three decades, the majority (i.e., over 70%) of the city was classified as Medium SUHI. Furthermore, the area of Very Low SUHI had moderate spatial changes and declined by about 2.26% from 1989 to 2019, which was mainly focused on the northern part of the study area. This may be mainly linked to the loss of vegetation covers (i.e., GS) and their replacement with urban impervious surfaces (i.e., BU). The area of Low SUHI increased by about 3.16% over the last three decades, of which further investigations revealed that 60.45% of this class was transited from the Very Low SUHI. Despite the variations of spatial distribution, the High SUHI maintained nearly 9% of Tehran's area over the study period. Moreover, Very High areas, except for some small patches, were mainly distributed and expanded in the western parts of Tehran, which was also in agreement with impervious surface developments. The Very High area increased by about 3.44% from 1989 to 2019 and reached 10.05% of the entire study area.   Table 3 for SUHI classes).
Additionally, the distribution of each LULC class at different SUHI categories for 1989, 1999, 2009, and 2019 was extracted by combining the LULC and classified SUHI maps. According to the results ( Figure 5), except for the WB, the majority area of other LULC classes had Medium SUHI. The results suggest that the SUHI intensity (i.e., SUHI category) has increased over the BL and BU areas over the last three decades while remaining relatively the same for the other two LULC classes. Moreover, the area of BU with High and Very High SUHI intensities slightly increased from 1989 to 2019, which indicates the necessity for precise monitoring of the SUHI phenomenon, especially in these regions, to adopt mitigation strategies. Likewise, one can see a notable increase in the area of BL with High and Very High SUHI intensities reaching nearly 47.53% in total over the last  Table 3 for SUHI classes).
Additionally, the distribution of each LULC class at different SUHI categories for 1989, 1999, 2009, and 2019 was extracted by combining the LULC and classified SUHI maps. According to the results ( Figure 5), except for the WB, the majority area of other LULC classes had Medium SUHI. The results suggest that the SUHI intensity (i.e., SUHI category) has increased over the BL and BU areas over the last three decades while remaining relatively the same for the other two LULC classes. Moreover, the area of BU with High and Very High SUHI intensities slightly increased from 1989 to 2019, which indicates the necessity for precise monitoring of the SUHI phenomenon, especially in these regions, to adopt mitigation strategies. Likewise, one can see a notable increase in the area of BL with High and Very High SUHI intensities reaching nearly 47.53% in total over the last three decades, which could also cause environmental issues for their proximities. The WB areas mostly included Very Low SUHI that was already expected based on the cooling effect characteristics of this land cover, except for a moderate fluctuation in 1999 that could be linked to the misclassification errors of the LULC maps. Although the area of Very High and High SUHI classes over the GS areas remained the same over time, the total area with Very Low and Low SUHI intensities declined by about 15.49% in 2019, compared to their corresponding utmost values (i.e., 39.12%) in 1999. In addition to the GS loss throughout the city, the thinning of GS cover can also be understood as another reason for this decline due to the weakening of cooling efficiency [76].
Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 23 throughout the city, the thinning of GS cover can also be understood as another reason for this decline due to the weakening of cooling efficiency [76]. To further appraise the impact of GS loss as a great source of SUHI mitigation [8,35], the decadal transition of GS to other LULC classes and the related SUHI intensity changes were investigated (see Table 4). Considering Table 4, the most significant transition of GS to other LULC classes, with a total area of 51.32 km 2 , occurred from 1999 to 2009. This considerable LULC change has also been accompanied by an increment in the SUHI intensity. In particular, the conversion of 20.61% of GS to BL from 1999 to 2009 resulted in a 2.16 °C rise in SUHI intensity in the corresponding locations. Likewise, the replacement of 18.61% of GS with BU areas led to a 1.24 °C growth in the SUHI intensity. After LULC changes between 1999 and 2009, the most extensive variation in GS occurred from 1989 to 1999 with a 40.03 km 2 loss of GS, which resulted in an average SUHI intensification of about 1.02 °C in the transition zones. In total, regarding the GS conversion to other LULC classes from 1989 to 2019, approximately 85.76 km 2 of GS disappeared, by which the SUHI intensity experienced an average rise of about 2.84 °C in the changed areas, leading to SUHI intensification throughout the city. To further appraise the impact of GS loss as a great source of SUHI mitigation [8,35], the decadal transition of GS to other LULC classes and the related SUHI intensity changes were investigated (see Table 4). Considering Table 4, the most significant transition of GS to other LULC classes, with a total area of 51.32 km 2 , occurred from 1999 to 2009. This considerable LULC change has also been accompanied by an increment in the SUHI intensity. In particular, the conversion of 20.61% of GS to BL from 1999 to 2009 resulted in a 2.16 • C rise in SUHI intensity in the corresponding locations. Likewise, the replacement of 18.61% of GS with BU areas led to a 1.24 • C growth in the SUHI intensity. After LULC changes between 1999 and 2009, the most extensive variation in GS occurred from 1989 to 1999 with a 40.03 km 2 loss of GS, which resulted in an average SUHI intensification of about 1.02 • C in the transition zones. In total, regarding the GS conversion to other LULC classes from 1989 to 2019, approximately 85.76 km 2 of GS disappeared, by which the SUHI intensity experienced an average rise of about 2.84 • C in the changed areas, leading to SUHI intensification throughout the city. Table 4. Mean Surface Urban Heat Island (SUHI) intensity of Green Space (GS) class along with the conversion rates to three other Land Use/Land Cover (LULC) classes based on the number of pixels, area per km 2 and the percentage of the converted pixels to the total number of GS pixels, and the mean SUHI intensity changes of converted pixels during the 10-year from 1989 to 1999, 1999

Relationship between UTFVI and LULC
In order to examine the status of the LULC classes in terms of environmental condition and thermal comfort, the relative percentage of each class for the two UTFVI classes of EC and WTC for 1989, 1999, and 2019 were extracted ( Figure 6). It is worth noting that only two extreme UTFVI classes were considered since these two classes cover nearly 90% of the study area. In general, during the last three decades, all LULC classes had downward trends in EC and upward trends in WTC. According to the results, BL areas were dominantly demarcated in WTC in all time intervals. The BL areas had the most environmental change in the last three decades. In particular, the area of BL with EC declined from 28.60% in 1989 to about 10% in 2019, and meanwhile, the BL area with WTC grew by about 22.53%, reaching 84.95% in 2019. In contrast, WB and GS areas were located mostly in EC; however, the GS experienced a more severe degradation (i.e., a downward trend in EC and an upward trend in WTC) in environmental conditions, which could probably be linked to the loss and thinning of GS cover. More importantly, Figure 6a shows that in 1989, over 61.35% of BU areas were located in EC, while in 2019, this percentage considerably decreased by about 19.01%, reaching 42.34%. Likewise, an unfavorable upward trend can be seen in Figure 6b, in which nearly 45.31% of the BU areas were demarcated in WTC. Consequently, the results confirmed that from 1989 to 2019, the BU, which is the permanent habitat of people, experienced significant thermal comfort decay. These adverse environmental conditions in BU areas can be associated with considerable urban expansions in unsuitable western parts of the city, which have also been classified as WTC in the earlier year of this study.
centage considerably decreased by about 19.01%, reaching 42.34%. Likewise, an unfavorable upward trend can be seen in Figure 6b, in which nearly 45.31% of the BU areas were demarcated in WTC. Consequently, the results confirmed that from 1989 to 2019, the BU, which is the permanent habitat of people, experienced significant thermal comfort decay. These adverse environmental conditions in BU areas can be associated with considerable urban expansions in unsuitable western parts of the city, which have also been classified as WTC in the earlier year of this study.

Intra-Annual Variation of SUHI
Seasonal variations of SUHI were also examined in the most recent year (i.e., 2019) of the study period to obtain an overview of the current spatio-temporal SUHI status. This analysis can provide practical information for future planning to promote a sustainable city. In this regard, the pattern and descriptive characteristics of the GSM are represented in Figure 7. Figure 7a-d illustrates the center and spatial footprint of primary SUHI, using an ellipse with an orientation of  and major and minor axes that are equal to a x and a y . As is clear, the SUHI throughout Tehran generally had higher intensity in summer and spring. It is evident that summer had the highest SUHI magnitude of 8.28 °C with the greatest SUHI footprint, while the winter had the lowest SUHI magnitude of 4.37 °C. The GSM results implied that almost the western and southwestern regions of the city, where

Intra-Annual Variation of SUHI
Seasonal variations of SUHI were also examined in the most recent year (i.e., 2019) of the study period to obtain an overview of the current spatio-temporal SUHI status. This analysis can provide practical information for future planning to promote a sustainable city. In this regard, the pattern and descriptive characteristics of the GSM are represented in Figure 7. Figure 7a-d illustrates the center and spatial footprint of primary SUHI, using an ellipse with an orientation of ϕ and major and minor axes that are equal to a x and a y . As is clear, the SUHI throughout Tehran generally had higher intensity in summer and spring. It is evident that summer had the highest SUHI magnitude of 8.28 • C with the greatest SUHI footprint, while the winter had the lowest SUHI magnitude of 4.37 • C. The GSM results implied that almost the western and southwestern regions of the city, where the center of SUHI was located in all seasons, had the highest rate of this phenomenon. This could be related to the development and concentration of industrial infrastructure and urban expansions in these regions. It is worth mentioning that there were other parts of the city with SUHI hotspots; however, the GSM model attempted to find and fit the greatest SUHI hotspot in the study area. Furthermore, the visual inspection of the results revealed that despite the concentrated patches with high SUHI intensity, the road networks with impervious surfaces, especially in the center of the city, included high SUHI intensity. This may be rooted in the fact that these roads in the city center mostly suffer from high traffic and the related heat from the combustion of fossil fuels.
As already discussed in detail in Section 4.1., the study area has undergone a significant LULC change in the past three decades. Here, the seasonal variations of SUHI intensities in different seasons between 1989 and 2019 were also explored. According to Table 5, the average SUHI intensities of all LULC classes in each season have modestly increased from 1989 to 2019. For instance, the BL average SUHI in summer 1989 increased by about 2.35 • C and reached 7.39 • C in summer 2019. Likewise, the average SUHI intensities over BU in spring from 1989 to 2019 rose from 2.86 • C to 4.44 • C. In total, the results revealed that the rise of the SUHI intensity itself through time and the LULC changes, in addition to causing a higher amount of LULC classes to experience higher SUHI intensities, have increased the average SUHI intensity throughout the city and decreased the thermal comfort and quality of the city.
of the city with SUHI hotspots; however, the GSM model attempted to find and fit the greatest SUHI hotspot in the study area. Furthermore, the visual inspection of the results revealed that despite the concentrated patches with high SUHI intensity, the road networks with impervious surfaces, especially in the center of the city, included high SUHI intensity. This may be rooted in the fact that these roads in the city center mostly suffer from high traffic and the related heat from the combustion of fossil fuels. As already discussed in detail in Section 4.1., the study area has undergone a significant LULC change in the past three decades. Here, the seasonal variations of SUHI intensities in different seasons between 1989 and 2019 were also explored. According to Table  5, the average SUHI intensities of all LULC classes in each season have modestly increased from 1989 to 2019. For instance, the BL average SUHI in summer 1989 increased by about 2.35 °C and reached 7.39 °C in summer 2019. Likewise, the average SUHI intensities over BU in spring from 1989 to 2019 rose from 2.86 °C to 4.44 °C. In total, the results revealed that the rise of the SUHI intensity itself through time and the LULC changes, in addition to causing a higher amount of LULC classes to experience higher SUHI intensities, have increased the average SUHI intensity throughout the city and decreased the thermal comfort and quality of the city.

Relationship between SUHI and AP
Air pollution is an influential factor on LST and consequently on SUHI intensity [77]. However, contradicting reports on the relationship between air pollutants and SUHI intensity in different regions [45][46][47][48] necessitate local studies to appraise their relationship in Tehran. In particular, these types of investigations should be examined through casespecific studies to ensure the selection of efficient strategies due to the complex and different characteristics of urban areas [46]. Accordingly, the results can assist in adopting more efficient strategies for the mitigation of these two adverse environmental phenomena [78]. Here, five different air pollutants, including NO 2 , O 3 , PM 2.5 , and PM 10 , were considered to discover their relationships with SUHI intensity. Figure 8 illustrates the generally positive relationships between daily pollutants concentration and SUHI intensity in 2019. To estimate the strength of bivariate associations between SUHI and air pollutants, the Correlation Coefficient (CC) values were calculated. Considering the annual investigations, O 3 and NO 2 had the highest and lowest CC values of 0.43 and 0.23, respectively. Furthermore, the CC values were also computed in seasonal time scales. According to the results, the highest CC values for O 3 (0.64) and PM 10 (0.60) were observed in autumn, while NO 2 (0. 43) had the highest CC value in spring. Additionally, PM 2.5 was the only pollutant that had almost steady behavior in all seasons, with an average CC value of 0.48. The discrepancies in seasonal CC values could be linked to other climate and meteorological factors such as wind speed, boundary layer structure, weather conditions, solar radiation, temperature inversion, and humidity [48,79]. The results indicate that these two environmental issues had a positive relationship, and an increase in one can strengthen the other and, thus, should be properly monitored since they are threats to human life and degrade the quality of life in metropolitan areas [80,81].

Discussion
Due to socio-economic and political conditions, Tehran has experienced a high rate of population rise in recent decades, which has led to Tehran being named one of the cities with the highest annual growth in Asia [34]. These demographic changes were accompanied by an increase in demands for urban facilities and resulted in a massive conversion of natural land covers to urban impervious surfaces [28]. This study illustrated that from 1989 to 2019, the dramatic urban expansion was accompanied by a considerable increase in the agglomeration of BU areas during each decadal interval (Figure 3). The sharp rise in BU constructions, especially in the western part of Tehran, has altered the thermodynamic properties of the surface and increased LST. The decadal examination of SUHI patterns throughout the city indicated the deterioration of thermal conditions in the study area. One of the main consequences of urbanization that emerged with the land cover transition was creating and intensifying SUHIs [82]. Some of the factors that have contributed to the aggregated growth of SUHI intensity include population growth, an increase in BU areas, the expansion of industrial activities, GS cover reduction, and changes in the intensity of solar radiation [52,83]. Different studies have focused on the appraisal of LULC transitions on increasing LST and SUHI. Furthermore, researchers have reported a negative relationship between GS cover and SUHI intensity and its significant role in mitigating the SUHI intensity [84]. Considering Tehran, the northern part of the city had the most considerable cooling effect associated with dense GS cover. However, over the last three decades, notable GS areas disappeared, and as a result, the SUHI intensity increased. Despite the GS loss, the thinning of GS density due to land cover changes also weakened

Discussion
Due to socio-economic and political conditions, Tehran has experienced a high rate of population rise in recent decades, which has led to Tehran being named one of the cities with the highest annual growth in Asia [34]. These demographic changes were accompanied by an increase in demands for urban facilities and resulted in a massive conversion of natural land covers to urban impervious surfaces [28]. This study illustrated that from 1989 to 2019, the dramatic urban expansion was accompanied by a considerable increase in the agglomeration of BU areas during each decadal interval (Figure 3). The sharp rise in BU constructions, especially in the western part of Tehran, has altered the thermodynamic properties of the surface and increased LST. The decadal examination of SUHI patterns throughout the city indicated the deterioration of thermal conditions in the study area. One of the main consequences of urbanization that emerged with the land cover transition was creating and intensifying SUHIs [82]. Some of the factors that have contributed to the aggregated growth of SUHI intensity include population growth, an increase in BU areas, the expansion of industrial activities, GS cover reduction, and changes in the intensity of solar radiation [52,83]. Different studies have focused on the appraisal of LULC transitions on increasing LST and SUHI. Furthermore, researchers have reported a negative relationship between GS cover and SUHI intensity and its significant role in mitigating the SUHI intensity [84]. Considering Tehran, the northern part of the city had the most considerable cooling effect associated with dense GS cover. However, over the last three decades, notable GS areas disappeared, and as a result, the SUHI intensity increased. Despite the GS loss, the thinning of GS density due to land cover changes also weakened the cooling efficiency of the remnant GS areas [76]. Taking the current structure of Tehran and the urban morphology, it may be practically infeasible and resource-intensive to create new large GS (i.e., parks and forests) areas within the hotspot zones of SUHI and, thus, it would be more cost-effective to establish new policies to encourage the use of green and cool roofs to mitigate the SUHI intensities [85], which can be followed up in future studies. Furthermore, utilizing high-albedo materials for roofs, pavements, and roads can also help mitigate SUHI and cool the city [86]. Additionally, the decadal analyses of the thermal comfort of the city through the UTFVI index demonstrated that the environmental condition of the city experienced a notable decay. More importantly, currently, 45.31% of BU exists in WTC zones that can potentially threaten human health. Despite other natural and climate-related factors that contributed to the environmental decay of Tehran, the construction of BU areas in locations with WTC in earlier decades has caused the unsuitable environmental condition of the majority of BU areas in Tehran. Therefore, it could be stated that the current status of the environmental condition of Tehran should be considered before any possible urban expansion to avoid further degradations of thermal comfort in BU areas.
Tehran is suffering from high-density air pollution [44]. Population growth, fossil fuel consumption, and the reduction of local winds due to vast vertical urban expansion have played essential roles in increasing the air pollutant concentration in Tehran. Additionally, urban growth, the conversion of natural land cover to buildings and human residents, and climate change have all enhanced the UHI intensity throughout the city. The preliminary analysis in this study reveals the positive interaction between SUHI and AP, which is also in accordance with similar studies in other areas [77,87,88]. In fact, high pollutant concentrations can trap more earth-emitted infrared radiation and heat in the urban environment, thus increasing the temperature [89]. Consequently, it is possible to incorporate suitable strategies to simultaneously reduce the SUHI intensity and AP concentration [90]. For instance, urban greening (i.e., increasing the vegetation cover through the city, especially in hotspot zones) is known to be an effective approach to mitigate the SUHI intensity. In this case, an efficient strategy would be choosing appropriate plant/tree species that have higher pollutant tolerance [90]. For example, Currie and Bass [91] found that extending green roofs by grass or shrubs can significantly mitigate AP, and it is well known that green roofs were stated to be one of the suitable approaches for UHI mitigation. Likewise, a model estimated by Yang et al. [92] indicated that increasing green roofs by about 20 ha can remove nearly 1700 kg of air pollutions during one year in Chicago. Furthermore, alterations in the urban microclimate can influence pollutant dispersion [90]. In particular, SUHI-related warm temperatures in urban areas affect the atmospheric turbulence and increase the boundary layer's height, leading to higher air pollution concentration [7,93]. Accordingly, SUHI mitigation will remove the effect that excessive warm temperatures have on atmospheric turbulence. For instance, Epstein et al. [94] reported that a decrease of about 0.35 k in the UHI intensity can decrease the urban boundary layer height by about 60%, thus boosting the dispersion of aerosol particles inside a deeper boundary layer, decreasing the aerosol concentration.
In this paper, Landsat-5 and Landsat-8 satellite images were collected from GEE. This geospatial cloud platform helps to reduce the dedicated time and effort for repetitive tasks, including pre-processing, image acquisition, and calibrating, of time-series satellite data [95]. Accordingly, it is recommended to exploit the considerable potential of GEE for long-term urban studies, especially for LULC, SUHI, and UTFVI mapping and monitoring.

Conclusions
Time-series Landsat-5 and Landsat-8 satellite images were employed to study the spatio-temporal variations of LULC, SUHI, and UTFVI throughout Tehran. The main findings of this study include the following: The generated LULC maps indicated that the simultaneous loss of GS areas along with the expansion of BU areas have contributed toward SUHI intensification by an average of 2.02 • C over the last three decades. The seasonal SUHI investigations revealed that the summer had the highest SUHI magnitude and footprint in Tehran. The UTFVI analysis revealed that the thermal comfort condition of Tehran has been significantly degraded since 52.35% of the city was identified as WTC, which was increased by about 16.37% compared with 2009. Notable thermal comfort degradation, especially in BU areas, was also associated with arbitrary urban expansion in the western part of Tehran since these regions were also classified as WTC even in earlier years when no BU existed. The preliminary results of comparing SUHI intensities and AP (i.e., NO2, O3, PM2.5, and PM10) concentrations indicated positive CC values between these two environmental phenomena in Tehran, which could be considered to adopt two-way mitigation strategies.
Policymakers and urban managers can employ the results of this study to avoid further environmental degradation in Tehran and assist with movement toward a sustainable city. In this regard, urban greening (i.e., by expanding urban parks, forests, and green roofs) can be considered to improve the environmental condition of Tehran. Determining and allocating pollution-tolerant plant species with reasonable cooling effects can be considered as the next step to enhance the quality of life in Tehran.    The NDVI images were calculated to apply emissivity correction of LST images. In total, 47 Landsat-5 and 21 Landsat-8 scenes were considered. There are several common image dates between the two above categories. Table A3. ImageCollection IDs of Landsat-5 and Landsat-8 datasets of this study.