Surface Urban Heat Island Assessment of a Cold Desert City: A Case Study over the Isfahan Metropolitan Area of Iran

Abstract

This study investigates the diurnal, seasonal, monthly and temporal variation of land surface temperature (LST) and surface urban heat island intensity (SUHII) over the Isfahan metropolitan area, Iran, during 2003–2019 using MODIS data. It also examines the driving factors of SUHII like cropland, built-up areas (BI), the urban–rural difference in enhanced vegetation index (ΔEVI), evapotranspiration (ΔET), and white sky albedo (ΔWSA). The results reveal the presence of urban cool islands during the daytime and urban heat islands at night. The maximum SUHII was observed at 22:30 p.m., while the minimum was at 10:30 a.m. The summer months (June to September) show higher SUHII compared to the winter months (February to May). The daytime SUHII demonstrates a robust positive correlation with cropland and ΔWSA, and a negative correlation with ΔET, ΔEVI, and BI. The nighttime SUHII displays a negative correlation with ΔET and ΔEVI.

1. Introduction

More than half the world’s population lives in urban environments [1], and it is predicted that by 2050 the population of the world’s cities will reach about 75% [2,3]. Two out of every three people are expected to live in urban areas by that time [4], so the 21st century will be the century of cities [5]. However, the rising urban population has brought many problems, such as stagnation and high concentrations of air pollutants [6,7], changes in rainfall patterns [8], land erosion [9], urban flooding [10], streams [11], the creation of new habitats [12], and urban heat island effects [13,14].

With increasing urbanization and changing natural landscapes to impervious surfaces, cities are more inclined towards the absorption of solar radiation [15], decrease in evapotranspiration [16], increase in runoff [17], augment surface friction [18], and release of anthropogenic heat [19]. Urbanization causes an increase in temperatures around urban areas and affects the urban thermal environment, thereby increasing the risk to human health. During heatwave conditions, when there is a lot of sunshine and a lack of wind to provide ventilation and disperse the warm air, urban environments can enhance the thermal condition and cause severe heat-related morbidity and mortality. The most well-known climate reform index of cities in environmental science is urban heat island (UHI) [20]. It forms when the urban environment’s temperature is higher than surrounding rural areas and has destructive effects on the ecosystem [21]. Urban heat islands were classified into three broad groups according to their characteristics in various layers of the urban atmosphere [22]: surface urban heat island (SUHI) [23], canopy layer urban heat island for micro scale analysis (CLUHI) [24], and boundary layer urban heat island for mesoscale analysis (BLUHI) [25]. Surface urban heat island (SUHI) is dependent on land surface temperature (LST) and has a strong correlation with the orientation of the surface relative to the sun as well as land use and land cover [26,27] which can be measured by radiometers onboard aircraft or satellites [28,29]. The great advantage of SUHI is that it can be conveniently measured across large spatial domains for a large number of cities [30]. This makes it possible to compare UHI between cities in large regions to explore the role of urbanization in affecting SUHI [31,32]. Canopy layer urban heat island (CLUHI) is typically quantified by the air temperature data of urban and rural stations [33,34], which are influenced by building geometry and the nature of pavement materials. The boundary layer urban heat island (BLUHI) is governed by the heated air from the upstream urban areas and the basic canopy layer where the canopy’s warm urban island occurs [24].

The long-term trends of SUHI and CLUHI were studied together in 272 cities in the mainland of China using satellite data and station-based air temperature data [35]. The studies reveal that the trend of the nighttime SUHI were strongly related to the CLUHI, whereas the relationship between the trends during daytime were relatively weak. The comparison of SUHI and CLUHI were studied recently for a better understanding of their vertical structure over 366 global cities within various background climates [36]. The findings of the study show that the annual mean SUHI is higher than CLUHI by 1.1 ± 1.9 °C during the daytime and 0.3 ± 1.5 °C at nighttime in equatorial, warm temperate, and snow climates. In contrast, in arid regions SUHI is lower than CLUHI by 0.8 °C. The relationship of SUHI and air temperature were investigated across 145 global large cities for the period 2003–2013, based on MODIS data [37]. The relationship was unstable in space, and the correlation differed significantly in macro-climatic conditions and was affected by the difference in vegetation between urban and rural areas.

SUHI has been widely studied in recent decades by analyzing the land surface temperature (LST) with available thermal remote sensing datasets of various spatial resolutions. The peak in SUHI intensity (urban LST minus rural LST) is generally observed at night [28,38]. The intensity of SUHI is weak during the daytime, but it becomes more pronounced after sunset. Various thermal remote sensing data are used to estimate the SUHI intensity, such as Landsat TM/ETM+/TIRS and Moderate Resolution Imaging Spectroradiometer (MODIS), which play an important role in SUHI research [14]. Generally, SUHI has been studied at a fine-scale using Landsat data because of its higher spatial resolution [39,40,41]. Meanwhile, MODIS LST data are used in regional and global SUHI studies due to their broad spatial coverage and high temporal resolution [31,42,43,44]. SUHI has been extensively studied over the different cities of the world situated over different climatic conditions. For example, [45] studied the seasonal SUHI effect over Siberian cities and found the strongest relationships between the UHI and population (log P). Additionally, [46] explored the indicators for quantifying SUHI over European cities using MODIS, and their study revealed that the temporal aspects and indicator selection were important factors in determining the SUHI intensity. Moreover, [47] studied the SUHI phenomenon over Ahmedabad, India, and found the impact of rural land cover dynamics in estimating the SUHI variation over semi-arid cities. Additionally, some studies also analyzed the underlying influencing factors of the SUHI, such as vegetation, albedo, evapotranspiration, meteorological condition, and anthropogenic heat emissions [48,49,50,51,52]. For example, [53] studied the SUHI phenomenon over 419 global cities and found that the intensity of SUHI was significantly negative with the vegetation during the daytime, while an insignificant relationship was observed during the nighttime. Also, [43] studied the SUHI footprint over 302 Chinese cities, and their findings indicated an increase in anthropogenic heat emission and decrease of vegetation activities, and they asserted that surface albedos should take lead responsibility for the expansion of SUHI footprint area.

However, little attention has been paid to SUHI studies over the cities of cold desert climates. The SUHI over desert cities shows contrasting behavior to that of the cities of temperate and humid zones. An inversion of SUHI during daytime is seen, with the city core appearing cooling than the surrounding suburban area during daytime. The cooling effect over the arid and desert cities has been found over cities around the world [54,55,56]. For example, the findings of [57] suggest that the nearby green areas play an important role in the control of atmospheric temperatures and in the promotion of citizens’ health and fog in the metropolis of Isfahan. In [58] found that SUHI-influenced areas mostly lie in the northern and southern parts of the city, where the vegetation cover is very low. A recent study estimated that the Isfahan metropolis area is cooler than the suburbs during the daytime, but at night it is around 2 K warmer than its surroundings [59]. With reference to the temporal and spatial actions of the SUHI metropolis of Isfahan, it seems that changes in the moisture, albedo, and composition of the atmosphere of the urban environment play a major role in the creation of SUHI. The impermeable surface expansion in the city of Isfahan was studied with Landsat data over the last three decades in the hot and cold months. The findings revealed that impermeable surfaces have increased by 2.8 times during 1985–2015 and observed a negative association with normalized difference vegetation index (NDVI) during hot months, while in the cold months, the urban region was colder than the surrounding rural area [60]. However, all the existing previous studies were focused on the diurnal analysis of SUHI and its relationship with vegetation using a short span of data, and they neglect the influence of other potential drivers of SUHI. Moreover, a detailed diurnal, seasonal, monthly, and temporal study is necessary along with the important driving factors to comprehend the SUHI phenomenon over Isfahan city, using a larger time series of data. In spite of the detailed SUHI investigation, the present study has some limitations, as it does not include the CLUHI study due to the non-availability of air temperature data. The combine study of SUHI and CLUHI can provide a better understanding of the vertical distribution of UHI and more in depth understanding of its dynamics, which are sometimes responsible for the formation of deadly heat wave phenomenon.

Therefore, this research is a comprehensive study of the SUHI phenomenon over the cold desert city of Isfahan, Iran. The main objective of the present research includes (i) revealing the annual, seasonal, monthly, and temporally spatial and temporal variation of LST from 2003 to 2019, (ii) estimating the annual, seasonal, monthly, and temporal variation of SUHI intensity during 2003–2019, and (iii) examining the relationship between the SUHI intensity and its associated influencing factors, including cropland, built-up area, vegetation, evapotranspiration, and albedo over the study period.

2. Study Area

The study was located in Iran in the cold desert city of Isfahan, with a latitude and longitude of 51°40′ E, 32°39′ N (Figure 1), and altitude of 1590 m [58]. The southern and western regions of Isfahan are surrounded by the Zagros mountains, while the northern and eastern regions are full of fertile plains [61]. Thus, this physiography creates a significant difference in terms of climate between the eastern and western regions of Isfahan [46]. According to Koppen–Geiger’s climate classification, Isfahan has a characteristic dry climate [62]. Dry weather and very low rainfall are the prominent features of this classification, where the minimum and maximum temperatures are −10.6 °C and 40.6 °C, respectively. The mean annual precipitation over the city is 116.9 mm [63].

According to national censuses, after industrial development and the establishment of different factories around it, the population of Isfahan rose from 255,000 in 1956 to 1,791,000 in 2010 and experienced remarkable urbanization. The urbanization ratio rose in 1956 and 2006, respectively, from 44.2 to 83.5 [64] and the proportion of farmland decreased, while the proportion of urban areas increased significantly during this period, primarily due to the decrease in the agricultural region of Isfahan [65].

3. Materials and Methods

3.1. Data Sources

The land surface temperature (LST) product obtained from MODIS sensor onboard Terra and Aqua satellite is used in the study. The MODIS sensor consists of 36 spectral bands ranging from 0.4 to 14.4 µm, providing data of different spatial resolutions (250 m, 500 m, and 1 km). The overpass time of the Terra and Aqua is between 10:30 a.m. and 13:30 p.m. local solar time during daytime, while during nighttime, it is between 22:30 p.m. and 01:30 a.m. local solar time.

First, the MODIS derived an 8-day LST product from Terra (MOD11A2) and Aqua (MYD11A2) satellite (https://modis.gsfc.nasa.gov/data/dataprod/mod11.php, accessed on 27 March 2021), available at 1 km × 1 km is used in this study. This LST product is derived using the two thermal infrared channels, band 31 (10.78–11.28 µm) and band 32 (11.77–12.27 µm) of the MODIS sensor. The 8-day LST product is estimated using a generalized split-window algorithm under clear-sky conditions [66], and the accuracy of the product is also checked with in-situ measurements, proving a bias less than 0.5 K [67]. We used the latest version (006) that removed all the cloud-contaminated pixels to enhance the accuracy. Moreover, the coefficients look up table (LUT) in this version is also updated, which considers a wide range of surface and atmospheric conditions, especially extending the upper boundary for (LST-Ts-air) in arid and semi-arid regions. Based on the quality control (QC) flag value, only pixels with high-quality LST were selected to eliminate the effect of retrieval algorithm and cloud cover. Data from 2003 to 2019 (17 years) from both Terra and Aqua satellite during both daytime and nighttime were used in this present study. The seasonal LST of Isfahan city is prepared separately for the summer (April, May, and June) and winter (December, January, and February) period. Besides, the monthly LST images were also prepared for the city to analyze the monthly variation.

LULC maps were extracted from the yearly MODIS land use land cover product, MCD12Q1 (https://lpdaac.usgs.gov/products/mcd12q1v006/, accessed on 27 March 2021), at a spatial resolution of 500 m from 2003 to 2019 and classified the Earth’s surface according to the International Geosphere-Biosphere Project (IGBP). The land cover classification scheme identifies a total of 17 classes, which include 11 classes of natural vegetation, 3 classes of mosaic use, and 3 classes of non-vegetated land. Based on the official urban limits of the city metropolitan areas studied, all pixels of the urban category inserted in the respective boundaries were defined as the urban core of the metropolitan area. The LULC maps were utilized for two main purposes: first, for urban and rural boundary area delineation and second, to extract the two important parameters (cropland land and built-up area), which were dominating in the study region for considering as SUHI drivers.

Besides, three more SUHI affective drivers were also analyzed in this present study: vegetation, evapotranspiration, and albedo. As for LST, vegetation and albedo data were extracted from the MODIS product. The enhanced vegetation index (EVI) was obtained from the MOD13A1 (16-day composite) at a spatial resolution of 500 m. The main algorithm of EVI generation is based on the constrained view angle and maximum value composite. In this approach, the number of observations with the highest vegetation value are compared and the observation with the smallest view angle, i.e., closest to nadir view, is chosen to represent the 16-day composite cycle. These methods inevitably result in spatial discontinuities because disparate days can always be chosen for adjacent pixels over the 16-day period. Thus, adjacent selected pixels may originate from different days, with different sun-pixel sensor viewing geometries and different atmospheric and residual cloud/smoke contamination. The white sky albedo (WSA) data are extracted from the daily MCD43A3 bi-hemispherical reflectance product at a spatial resolution of 500 m. The MCD43A3 albedo model parameters product supplies the weighting parameters associated with the Ross Thick Li Sparse Reciprocal BRDF model that best describe the anisotropy of each pixel. These parameters can be used in a forward version of the model to reconstruct the surface anisotropy effects and thus correct directional reflectance’s to a common view geometry or to compute the integrated black-sky and white-sky albedos. Alternatively, these parameters can be used with a simple polynomial to easily estimate the black-sky albedo with good accuracy for any desired solar zenith angle. Data with the best (QC = 0) and good (QC = 1) quality flags were used. The evapotranspiration (ET) data used in this study are model data obtained from the operational simplified surface energy balance (SSEB) model [68]. This model provides monthly ET data at a spatial resolution of 1 km. The SSEB model has a unique parameterization for operational applications with the pre-defined, seasonally dynamic boundary condition that is unique to each pixel. The original formulation of this SSEB model is based on the hot/cold pixel principles of the SEBAL [69] and METRIC [70] models. SEBAL (Surface Energy Balance Algorithm for Land) uses surface temperature, hemispherical surface reflectance, and Normalized Difference Vegetation Index (NDVI), as well as their interrelationships to infer surface fluxes for a wide spectrum of land types. METRIC (Mapping Evapotranspiration at High Resolution using Internalized Calibration) uses the instantaneous to daily ET extrapolation method called alfalfa reference ET fraction, which employs wind speed and air temperature that according to better incorporate local/regional surface/environmental conditions than the evaporative fraction of other remote sensing ET algorithms.

3.2. Methods

3.2.1. Definitions of Urban and Rural Regions

The MODIS yearly land use land cover (LULC) data product (MCD12Q1) at a spatial resolution of 500 m is used in this study to delineate the urban and rural regions from 2003 to 2019. Firstly, a square box around the city of Isfahan is chosen in such a way that the urban and surrounding non-urban pixels are nearly in an equal ratio. The MODIS LULC data consist of five different classes like open shrubland, grassland, cropland, urban or built-up land, and barren land. Urban or built-up land is used to delineate the urban pixels, while all other four classes were used to delineate the non-urban pixels [71], as shown in Figure 2. As urban area changes continuously, in order to avoid biases in our result by using the same boundary, we changed our urban and rural boundaries in the years 2003, 2007, 2011 and 2015. The mean LULC area for the above years is shown in

Table A1. Mean area (km²) of the different LULC classes.

LULC	2003	2007	2011	2015
Open Shrub land	218.75	243.58	285.73	280.77
Grassland	76.27	75.10	80.81	102.34
Cropland	414.81	404.67	344.75	316.01
Urban or Built-up	319.86	326.46	339.53	350.55
Barren land	480.20	467.09	479.08	490.28

of Appendix A. Whereas, the temporal LULC maps from 2003 to 2019 is shown in Figure A1 of Appendix A.

3.2.2. Surface Urban Heat Island Assessment

The urban and non-urban (rural) boundary extracted from MCD12Q1 is used to calculate surface urban heat island (SUHI) intensity. The Simplified Urban-Extent (SUE) algorithm is used in this study to calculate to intensity of SUHI [71,72]. As per this algorithm, SUHI defines as the difference in the LST of the urban pixels and the non-urban pixels (Equation (1)).

$$\mathrm{SUHI}={\mathrm{meanLST}}_{\mathrm{urban}\mathrm{pixels}}-{\mathrm{mean}\mathrm{LST}}_{\mathrm{non}-\mathrm{urban}\mathrm{pixels}}$$

(1)

Figure 2 shows the steps for estimating the SUHI using 500 m MCD12Q1 LULC from a 1 km MODIS LST product. Firstly, the urban and non-urban extent are clipped from land cover data for the years 2003, 2007, 2011, and 2015 as defined in Section 3.2.1. After sub setting, the mean LST of both subsets are calculated for each year, and their difference is the surface UHI (SUHI) for that particular year.

3.2.3. LST Trend Analysis Method

We used the Mann–Kendall test to detect the temporal trend in the LST datasets [73]. The test has been extensively used in the hydro-climatic time series data and has always proven to be an efficient tool in comparison to other available tools for detecting trends [74,75,76,77,78]. This test has numerous advantages as it analyzes time series that are not required to follow a specific linear or nonlinear trend. We calculated the standardized Mann–Kendall statistics, which display whether there is a significant trend is present in the datasets or not at a specific significance level, p. A significance level of p = 0.5 is considered in this study. Besides, the magnitude of the trend is also determined by Sen’s slope estimator test [79]. A positive value indicates an increasing trend and a negative value indicates a falling trend.

3.2.4. Exploring the SUHI Drivers

It is very well known that the changes in land use land cover have significantly affected the seasonal and diurnal variation of SUHI [14]. In the present study, five different potential driving variables of SUHI are considered to understand the reason for SUHI variability over the study period. Two of them, cropland and built-up index, are extracted from MODIS LULC product, which predominated in the study area. The vegetation, evapotranspiration, and albedo are the other three important factors (extracted from MODIS) that can affect SUHI. The differences between an urban and rural area in enhanced vegetation idea (ΔEVI), evapotranspiration (ΔET), and white sky albedo (ΔWSA) were considered to investigate the SUHI variability. This difference was estimated exactly the same way as we did for SUHI estimation in Equation (1). ΔEVI represents the difference in mean EVI over rural pixels and mean EVI over non-rural pixels, and ΔET and ΔWSA are estimated in the same manner. Then, a Pearson’s correlation coefficient between the potential driving variables and SUHII across the city and the years was computed.

4. Results

4.1. Spatial Variation and Temporal Trend of LST

The diurnal, seasonal, annual, monthly, and temporal variation of LST from 2003 to 2019 are analyzed in this section. The daytime (10:30 a.m. and 13:30 p.m.) and nighttime (22:30 p.m. and 01:30 a.m.) variation was shown in this study.

4.1.1. Annual Spatial Pattern of LST

The spatial variation (first row) and temporal trend (second row) of the mean annual LST during both daytime and nighttime from 2003 to 2019 are shown in Figure 3. It is evident that the SUHI effect over the city was more pronounced during the nighttime compared to the daytime. The mean LST over urban (U) and non-urban (NU) pixels (as define in Section 3.2.1) was also calculated as in

Table 1. Annual, seasonal (summer and winter), and monthly (January-December) mean LST from 2003 to 2019 for urban and non-urban pixels at daytime (10:30 a.m. and 13:30 p.m.) and nighttime (22:30 p.m. and 01:30 a.m.).

		Terra 10:30 a.m.	Aqua 13:30 p.m.	Terra 22:30 p.m.	Aqua 01:30 a.m.
Annual	Non-Urban Pixels	29.9 ± 2.3	33.4 ± 2.5	10.1 ± 1.3	8 ± 1.4
	Urban Pixels	28.9 ± 1.9	32.6 ± 2	11.4 ± 1.7	9.2 ± 1.9
Summer	Non-Urban Pixels	37.4 ± 3	40.5 ± 3.2	15.4 ± 1.4	12.6 ± 1.6
	Urban Pixels	36.6 ± 2.3	40.2 ± 2.5	16.6 ± 1.7	13.8 ± 1.9
Winter	Non-Urban Pixels	14.5 ± 1.6	18.5 ± 1.8	−0.5 ± 1.1	−1.8 ± 1.1
	Urban Pixels	13.6 ± 1.3	17.7 ± 1.5	0.7 ± 1.6	−0.5 ± 1.8
Jane	Non-Urban Pixels	11.9 ± 1.6	16.1 ± 1.8	−1.9 ± 1.2	−3 ± 1.2
	Urban Pixels	11.1 ± 1.2	15.4 ± 1.4	−0.5 ± 1.7	−1.6 ± 1.8
February	Non-Urban Pixels	18.4 ± 1.9	23.1 ± 1.9	1 ± 1.1	−0.7 ± 1.1
	Urban Pixels	17.5 ± 1.6	22.2 ± 1.6	2.2 ± 1.6	0.5 ± 1.7
March	Non-Urban Pixels	25.9 ± 2	29.6 ± 2.1	5.7 ± 1.1	3.4 ± 1.2
	Urban Pixels	25.3 ± 1.5	29 ± 1.6	6.9 ± 1.6	4.5 ± 1.7
April	Non-Urban Pixels	30.1 ± 2.4	33.1 ± 2.5	9.9 ± 1.2	7.4 ± 1.3
	Urban Pixels	29.6 ± 1.8	32.9 ± 1.9	11 ± 1.5	8.4 ± 1.7
May	Non-Urban Pixels	36.3 ± 3.2	39.2 ± 3.3	14.7 ± 1.4	12 ± 1.5
	Urban Pixels	35.7 ± 2.4	39.1 ± 2.6	15.9 ± 1.6	13.2 ± 1.8
June	Non-Urban Pixels	43.9 ± 3.5	47.3 ± 3.6	20.1 ± 1.7	17 ± 1.9
	Urban Pixels	42.8 ± 2.8	46.6 ± 2.9	21.5 ± 1.8	18.3 ± 2.1
July	Non-Urban Pixels	45.8 ± 3.5	49.6 ± 3.7	22.2 ± 1.7	19 ± 1.9
	Urban Pixels	44.4 ± 3	48.5 ± 3.1	23.6 ± 1.8	20.2 ± 2.2
August	Non-Urban Pixels	43.6 ± 3.6	47.5 ± 3.9	20.1 ± 1.9	17 ± 2.1
	Urban Pixels	42.3 ± 3	46.4 ± 3.2	21.4 ± 2	18.2 ± 2.3
September	Non-Urban Pixels	38.4 ± 3.2	41.8 ± 3.4	15.9 ± 1.9	13.2 ± 2.2
	Urban Pixels	37 ± 2.7	40.6 ± 2.8	17.2 ± 2.1	14.4 ± 2.3
October	Non-Urban Pixels	30.7 ± 2.5	33.3 ± 2.5	10.3 ± 1.7	8 ± 1.9
	Urban Pixels	29.3 ± 2.2	32.1 ± 2.2	11.4 ± 1.9	9 ± 2.1
November	Non-Urban Pixels	20.3 ± 1.7	22.7 ± 1.9	3.2 ± 1.3	1.8 ± 1.2
	Urban Pixels	19.2 ± 1.6	21.4 ± 1.7	4.3 ± 1.7	2.7 ± 1.7
December	Non-Urban Pixels	13.1 ± 1.5	16.2 ± 1.7	−0.8 ± 1.2	−1.7 ± 1.2
	Urban Pixels	12.3 ± 1.2	15.3 ± 1.5	0.5 ± 1.7	−0.4 ± 1.8

. During the daytime, the mean annual LST was higher in non-urban pixels compared to urban pixels at both 10:30 a.m. and 13:30 p.m. The observed LST is higher at 13:30 p.m. (33.4 ± 2.5 °C for NU and 32.6 ± 2 °C for U pixels), than 10:30 a.m. (29.9 ± 2.3 °C for NU and 28.9 ± 1.9 °C for U pixels). In contrast, during nighttime, a reverse phenomenon is observed, witnessing a higher LST over urban pixels compared to non-urban pixels at both 22:30 p.m. and 01:30 a.m. The mean LST was lower at 01:30 a.m. (8 ± 1.4 °C for NU and 9.2 ± 1.9 °C for U pixels) in comparison to 22:30 p.m. (10.1 ± 1.3 °C for NU and 11.4 ± 1.7 °C for U pixels). The surrounding land cover is dominated by barren land, which can absorb heat much more rapidly than built-up land after sunrise, explaining the higher LST in non-urban pixels. Meanwhile, urban built-up land can retain heat for longer time than barren land cover due to the higher thermal inertia, resulting in the slow emission of heat from built-up lands. Due to this reason, the urban pixels witness a higher LST than the non-urban pixels in nighttime.

The second row in Figure 3 shows Sen’s slope trend of mean annual LST during 2003–2019. During the daytime, the central urban area has a lower trend in mean LST compared to the surrounding rural area. However, a significantly dispersed increasing trend is observed during the nighttime compared to the cluster trend in the daytime. The magnitude of Sen’s slope was higher in the daytime (0.49 to −0.26 °C/year) compared to the nighttime (0.14 to −0.05 °C/year) trend.

4.1.2. Seasonal Spatial Pattern of LST

The seasonal variation of mean LST during 2003–2019 is shown in Figure 4a,b, where the first row represents the mean LST, and the second-row represents Sen’s slope trend during the study period. The daytime represents a clear heterogeneity in the LST distribution compared to the nighttime imagery during both the season. The summer season shows higher land surface temperature (50.2 to 31.2 °C during daytime and 19.5 to 9.2 °C during the night) compared to the winter season (24.1 to 7.3 °C during daytime and 3.6 to −6.9 °C during the night) as shown in Figure 4a,b respectively. More clear evidence of the surface urban heat island phenomenon is visible in the nighttime, while in the daytime, a cool island effect is seen, where the surrounding non-urban regions show a higher LST than the central urban region. The mean LST over urban and non-urban pixels were also calculated as shown in Table 1. During daytime, the mean LST over urban pixels (36.6 ± 2.3 °C at 10:30 a.m. and 40.2 ± 2.5 °C at 13:30 p.m. during summer and 13.6 ± 1.3 °C at 10:30 a.m. and 17.7 ± 1.5 °C at 13:30 p.m.) is comparatively lower than that of non-urban pixels (37.4 ± 3 °C at 10:30 a.m. and 40.5 ± 3.2 °C at 13:30 p.m. during summer and 14.5 ± 1.6 °C at 10:30 a.m. and 18.5 ± 1.8 °C at 13:30 p.m.). Meanwhile, at night, the urban pixels show a higher LST (16.6 ± 1.7 °C in summer and 0.7 ± 1.6 °C in winter at 22:30 p.m., and 13.8 ± 1.9 °C in summer and −0.5 ± 1.8 °C in winter at 1:30 a.m.) compared to non-urban pixels (15.4 ± 1.4 °C in summer and −0.5 ± 1.1 °C in winter and 22:30 p.m., and 12.6 ± 1.6 °C in summer and −1.8 ± 1.1 °C in winter at 1:30 a.m.), witnessing the presence of the urban heat island phenomenon.

The result of Sen’s slope suggests a higher increasing trend in the summer season compared to the winter season which is also higher in the daytime than the nighttime within the same season. The center urban area witnesses a lower trend compared to the surrounding rural area in the summer season (second row of Figure 4a), while the winter season shows a relatively different trend pattern than summer and annual. In winter nighttime, a positive LST trend is observed in an urban area (more pronounced at 1:30 a.m.) as compared to daytime (Figure 4b). The trend varies from 0.57 to −0.35 °C/year in the daytime to 0.25 to −0.04 °C/year in the nighttime in the summer season. By contrast, in the winter season, it varies from 0.37 to −0.19 °C/year in the daytime to 0.23 to −0.10 °C/year in the nighttime.

4.1.3. Monthly Spatial Pattern of LST

The monthly spatial variation of mean LST during 2003–2019 over the city is also analyzed in this study during both daytime and nighttime (Figure 5). The result obtained from monthly variation also suggests that the SUHI was intensified during the nighttime (more at 22:30 p.m. as compared to 1:30 a.m.) compared to during daytime, in which an urban cool island effect can be seen. The mean LST variation is plotted as a box and whisker plot in Figure 6. It can be seen that the maximum mean LST was observed during the month of July, while the minimum was observed in January. Noticeably, the mean LST at 13:30 p.m. was always higher than at 10:30 a.m. during the daytime, while during nighttime, it was always lower at 1:30 a.m. compared to 22:30 p.m. This also suggests an occurrence of the urban cool island during daytime as non-urban pixels had a higher mean LST than urban pixels at 10:30 a.m. and 13:30 p.m. However, during the nighttime (both 22:30 p.m. and 1:30 a.m.), the mean LST over urban pixels was always higher than non-urban pixels suggesting a positive urban heat island over the city (Table 1).

Further, the pixel-based Sen’s slope trend of mean LST is also calculated at a monthly level (Figure 7) from 2003 to 2019 during both the daytime and nighttime. The trend was comparatively higher in the daytime hours compared to nighttime. During the daytime, the maximum trend was observed from June to September (in the range of 0.70 to 0.87 °C/year), while the minimum trend was observed in the months of November to February (in the range of 0.19 to 0.36 °C/year. At nighttime, the maximum trend was only 0.33 °C/year in September and the least in February with 0.16 °C/year.

4.2. SUHI Intensity Variation

The annual, seasonal, and monthly variations of SUHI intensity during daytime and nighttime hours are shown in Figure 8. The results clearly indicate the presence of urban cool islands during the daytime and urban heat islands in the nighttime. During the daytime, the SUHI intensity was much lower at 10:30 a.m. compared to at 13:30 p.m. The months from July to October show minimum SUHI intensity values compared to other months. In the nighttime, there was not much variation in SUHI values thought out the year. The intensity of SUHI was higher at 22:30 p.m. compared to that at 1:30 a.m.

4.3. Relationships between the SUHI Intensity and Its Potential Influencing Factors

The relationship between SUHI intensity (SUHII) and its potential influencing factors was examined, which contribute to understanding the specific cause of the SUHI effect. To explore the cause of SUHI variability, five different variables were analyzed and correlated with SUHII: Cropland area (Crop), built-up area (BI), ΔET, ΔEVI, and ΔWSA (i.e., ET, EVI, and WSA urban–rural difference). Figure 9a shows the annual, seasonal (summer and winter), and monthly (January to December) Pearson’s correlation coefficient plot of SUHII with the different driving parameters, average to yearly level during daytime (Terra 10:30 a.m. and Aqua 13:30 p.m.), while Figure 9b shows during nighttime (Terra 22:30 p.m. and Aqua 01:30 a.m.). The daytime SUHII shows a significant strong positive correlation with Crop and ΔWSA, while negative with ΔET, ΔEVI, and to some extent BI as well. The daytime correlation was more significant at 13:30 p.m. than 10:30 a.m. The nighttime SUHII shows very less correlation except with the negative correlation of ΔET and ΔEVI. For seasonal correlation, the summer has a much more significant correlation with the driving parameters compared to the winter.

The results of monthly variation reveal that the daytime SUHII was significantly positively correlated with Crop irrespective of the months. The BI shows a significant correlation with SUHII during daytime (more so in the months of July, August, and September), while in the nighttime, the correlation was very less and insignificant. ΔET shows a significant negative correlation in the daytime compared to the nighttime, especially in the months from March to September. The nighttime correlation of ΔET with SUHII is less negative than daytime, where a slight positive correlation is seen from January to March. The correlation of ΔEVI was also negative with SUHII and more in the daytime than nighttime, and higher in the summer months (March–May) than in the winter months. The correlation of ΔWSA is the opposite, as it shows a positive correlation during daytime (higher in summer months than winter months) and negative in the nighttime (significant in winter months, especially in November and December).

5. Discussion

5.1. Variations in the SUHI Intensity

The present study examines the diurnal, annual, seasonal, and monthly behavior of LST and SUHI intensity over Isfahan city from 2003 to 2019 using MODIS datasets. The diurnal and monthly variation of the SUHI index can be found in the previous study [59]. Their findings suggest that the Isfahan metropolitan area is 3–4 K colder than the surrounding suburban area during the daytime, while at nighttime, it is 2 K warmer than its surroundings. The result of the present study also found evidence of the occurrence of the urban cool island during daytime 10:30 a.m. and 13:30 p.m.), while at nighttime (22:30 p.m. and 01:30 a.m.), a positive heat island phenomenon is profoundly visible over the city. The SUHI intensity was found to be minimum at 10:30 a.m. and maximum at 22:30 p.m. The lower land surface temperature over the Isfahan city as compared to surrounding rural area during daytime is also found in previous studies using Landsat data [57,58,60,80]. Moreover, another study [57] suggests that the central portion of the Isfahan metropolis exhibits the highest portion of cold and very cold classes of LST. The findings also suggest that the northern, southern, and eastern portions of the city show much warmer LST, which is also seen in the present study. Considering the seasonal behavior of SUHI, it is evident that the SUHI intensity is weaker in the summer season and stronger in the winter season [59]. The maximum SUHI intensity was found during the January month while the minimum in the months of the summer season.

The phenomenon of a higher LST over the urban area compared to the surrounding rural areas pertaining to negative SUHI during daytime is a typical characteristic of cities situated over arid regions [14,54,55,81]. Moreover, during the daytime, there was an obvious seasonal variation in SUHI, except for the city of the equatorial zone. The amount of rural vegetation cover, direct solar radiation, and longer sunshine duration in summer may contribute to the seasonal variation in SUHI [82]. A negative SUHI during daytime is also reported in various urban areas of the world, including northwestern China [83], central Asia [53], western United States [84], and cities in India [85,86,87] and South America [88]. All these urban areas are situated over arid or semi-arid regions, a similar climatic characteristic prevailing over the present study area of Isfahan city. Cities situated over arid and semi-arid climatic zones such as Abu Dhabi, Kuwait, Riyadh, Las Vegas, and Phoenix experienced lower LST than the suburban area during daytime. The data suggest that this negative SUHI is due to the existence of a large number of bare lands present in suburban areas, which mostly absorb sunlight rather than reflecting it, leading to a higher temperature in suburban areas than in city cores [59,89,90]. The present study area’s dominating land cover is bare land in the suburban region, which leads to a higher land surface temperature than the urban area of the city.

5.2. The Effects of Each Factor on SUHI Intensity

The second objective of the present research was to explain the potential driving variables affecting the diurnal, seasonal, and monthly variation of SUHI intensity. The SUHI phenomenon is driven by various factors that can be broadly categorized as LULC distribution and its spatial pattern, urban site characteristics, and landscape configuration. During the daytime, SUHII was negatively correlated with ΔET, meaning that an increase of SUHI intensity is related to a reduction of the urban ET or a rise of the rural ET. The evaporation from the water surfaces and transpiration from green leaves provides cooling to the rural area, resulting in positive SUHI over the cities. The evaporative cooling is a much more dominant controlling parameter in the daytime for cities over dry climatic regions [55]. The results of the study conducted by [86] suggest that the ET increase in the city area due to higher water use and gardening with irrigating resulted in negative SUHII during the daytime. The relationship between the SUHI intensity and vegetation was significant during the daytime due to the evaporative cooling affect by vegetation [53,55]. The correlation for summer days was significant and negative between the SUHII and ΔEVI in this study compared to winter days, which was similar to the results of the other previous studies [47,53,91]. The vegetation over the city is less affected by seasons, and so a decreasing ∆EVI can intensify the SUHII, even in the summer. The lack of vegetation transpiration in the absence of sunlight during the nighttime leads to an insignificant relation of SUHII with ΔEVI [91,92]. The rural area with higher vegetation cover can increase the latent heat flux via the process of transpiration, which provides cooling, thus lowering the LST of the urban area and resulting in a cool urban heat island during the daytime. It is inferred from the previous studies that vegetation plays an important role in the mitigation of the SUHI effect [93]. The study in Guangzhou, China, found that a decrease in 16% of urban vegetation from 1990–2007 caused an increase in LST by 2.5 °C [94]. Except for ΔEVI and ΔET, there is no significant negative correlation of SUHI intensity with other parameters. Cropland and ΔWSA are found to be effective SUHI controlling variables during the nighttime with a significant positive correlation. In the winter, the LST over cropland is lower than the urban LST due to irrigation and plantation growth, and the rural LST is lower during nighttime, which causes an insignificant correlation with SUHI compared to the summer season. Urbanization generally decreases surface albedo and emissivity, largely due to the urban canopy effect and the nature of the newly added urban material [95,96]. They are responsible for reducing heat loss through increased heat storage and net all-wave radiation in urban areas. The night SUHII is mainly driven by the surface heat fluxes that originate from heat storage during daytime, resulting in the positive correlation of SUHII with the urban–rural difference in albedo (ΔWSA) in the nighttime as seen in our study. This result is also supported by the previous research that suggests urban regions with lower surface albedo tend to increase the diurnal LST and have more energy for releasing at night [97]. Moreover, the recent study on SUHI study also supports a strong diurnal control of albedo in SUHI variation [55].

6. Conclusions

The present study investigated the diurnal, seasonal, annual, and monthly spatiotemporal variation of the LST and SUHII of Isfahan city, Iran, during the last 17 years (2003–2019) based on the MODIS datasets of 1 km resolution. Correlation analyses were also conducted to reveal the relationship between the SUHII and its driving factors, including cropland, built-up area, ΔEVI, ΔET, and ΔWSA.

The results show the presence of an urban cool island during the daytime and a positive heat island during the nighttime. Considering the diurnal variation of SUHI, the maximum SUHII was observed at 22:30 p.m., while the minimum was observed at 10:30 a.m. The SUHII values ranged from −1.5 to +1.5 °C over the metropolis area of Isfahan. This seasonal variation suggests that the SUHI was much more dominating in the winter season compared to the summer season, with both higher and lower SUHII value during daytime and nighttime in the winter season. The SUHII were higher from June to August months, while minimum from February to May months. The SUHI intensity was significantly and negatively correlated with ΔEVI and ΔET during the daytime. This negative correlation suggests that the negative SUHI during daytime is ascribed to the low vegetation activity in the rural area, dominated by croplands and bare land, which absorb much more heat and warm the surrounding rural surfaces. Furthermore, the nighttime SUHII was related to cropland and ΔWSA. The lower albedo material in the urban area due to urbanization is responsible for reducing heat loss through increased heat storage and net all-wave radiation in urban areas. The nighttime SUHI mainly drivees the surface heat flux driven by heat storage. Therefore, a strong positive correlation of WSA is observed with the SUHII over Isfahan during the nighttime.

Due to the significant SUHI effects observed over Isfahan city in recent years, variation in the urban thermal environment and its related ecological responses should be comprehensively analyzed. The result of present study can be an important reference for understanding the spatiotemporal variations in the SUHI effect, and the interaction between human activities and land-surface ecosystems. However, some uncertainties remain, and other factors (e.g., soil moisture, thermal inertia, landscape configuration) that are associated with the SUHI phenomenon should be thoroughly studied in the future work.