Integrated Influencing Mechanism of Potential Drivers on Seasonal Variability of LST in Kolkata Municipal Corporation, India

Abstract

Increasing land surface temperature (LST) is one of the major anthropogenic issues and is significantly threatening the urban areas of the world. Therefore, this study was designed to examine the spatial variations and patterns of LST during the different seasons in relation to influencing factors in Kolkata Municipality Corporation (KMC), a city of India. The spatial distribution of LST was analyzed regarding the different surface types and used 25 influencing factors from 6 categories of variables to explain the variability of LST during the different seasons. All-subset regression and hierarchical partitioning analyses were used to estimate the explanatory potential and independent effects of influencing factors. The results show that high and low LST corresponded to the artificial lands and bodies of water for all seasons. In the individual category regression model, surface properties gave the highest explanatory rate for all seasons. The explanatory rates and the combination of influencing factors with their independent effects on the LST were changed for the different seasons. The explanatory rates of integration of all influencing factors were 89.4%, 81.4%, and 88.7% in the summer, transition, and winter season, respectively. With the decreasing of LST (summer to transition, then to winter) more influencing factors were required to explain the LST. In the integrated regression model, surface properties were the most important factor in summer and winter, and landscape configuration was the most important factor in the transition season. LST is not the result of single categories of influencing factors. Along with the effects of surface properties, socio-economic parameters, landscape compositions and configurations, topographic parameters and pollutant parameters mostly explained the variability of LST in the transition (11.22%) and summer season (15.22%), respectively. These findings can help to take management strategies to reduce urban LST based on local planning.

1. Introduction

World urbanization is increasing, with more than 55% of the world’s population living in cities, and this figure is expected to reach up to 68% by 2050 [1]. Due to rapid urbanization, the world’s urban area has increased by 168% from 2001 to 2018, with the largest growth in Asia and Africa [2]. Over the past three decades, India has experienced rapid urban growth due to the new economic reform during the 1990s [3]. In 2020, more than 35% of the population was residing in urban areas [4], but almost 50% of India’s population will live in urban areas by 2050 [5]. The rapid urbanization process will transform the natural landscape into an impervious surface, which will alter the radiative, thermal, moisture, and emission properties [6], thereby leading to a significant increase in land surface temperature (LST) [7]. Due to the increase in urban LST, extreme weather events, such as the frequency and intensity of heat waves and extreme summer rainfall, have increased in urban areas [8]. Furthermore, the energy consumption rate [9], air pollution [10], health risks, and heat-related deaths [11,12] have increased throughout the world’s urban areas. Therefore, there has been increasing interest among researchers of urban LST for the development of suitable strategies for sustainable development and to create an eco-friendly environment in urban areas [13].

Generally, three major types of temperature are considered: Urban Canopy Layer (UCL), Urban Boundary Layer (UBL), and Urban Surface Layer (USL) [14]. UCL represents the thermal environment from the surface to mean building height [14], and measurements are based on the data from meteorological stations, field measurements, and model simulations [15]. UBL represents the thermal environment above the UCL [16], generally based on special platforms, such as radiosondes, aircraft, weather stations, and flux towers [14]. USL represents the Land Surface Temperature (LST), which is primarily based on in situ measurements from ground stations [17] and remote sensing-based measurements from satellite images providing a more representative study [18]. Although in situ-based measurements provide higher spatial resolution with better accuracy, but the inadequate distribution of stations and small area coverage complicates their use [19]. To overcome the inadequate small area coverage and to cover the large areas, satellite-based studies are increasing the attention around estimating urban LST [20]. Moreover, satellite images provide various degrees of spatial resolution with large-scale spatial coverage, which promotes the use of remote sensing-based studies [21].

In heterogeneous urban environments, spatial patterns and trends of LST depend on various influencing factors [22], and the factors vary from one city to another [23]—as well as varying seasonally [24]. Thus, it is important to investigate the factors of urban LST during the different seasons and how LST is quantitatively affected by influencing factors, and which factors are very important. The various previous studies influencing factors of LST are divided into certain categories, namely, landscape configurations [25,26], socio-economic indicators [27], biophysical variables [28], and landscape compositions [24,26]. However, different topographic parameters, such as slope, elevation, aspect, and hill shade, affect LST [29]. Furthermore, air pollutants, such as aerosol optical depth [30,31] and greenhouse gasses concentration (NO₂, SO₂, CO), could also affect urban LST [30,32]. Therefore, a numerous amount and categories of influencing factors affect LST, but most of the studies analyzed one or only a few categories of influencing factors in their study. Previous studies [7,22,24,25,26,27,28] did not incorporate all these categories of influencing factors to estimate the integrated effects of these factors on LST during different seasons. Furthermore, prior studies failed to analyze the effects of each category separately and to identify the dominant factors of each category of influencing factors during the different seasons.

To address these problems, Kolkata Municipality Corporation (KMC) (for details see Section 2.1 Descriptions of study area) was used as the study area. Various previous studies have been conducted on the patterns and trends with relation to few limited influencing factors in KMC areas [33,34]. However, KMC consists of 141 wards, and their surface cover is different in each ward. Each surface cover has properties that have unique radiative, thermal, moisture, and aerodynamic properties [35]. Therefore, it is important to study each ward’s spatial variation of LST. However, previous studies have not addressed the spatial variation of LST in the different wards. Furthermore, little is known regarding the variation of LST in the different LULC types and influencing factors for the different seasons.

For these consequences, this study was designed to explore the spatial variations in LST with relation to the influencing factors during the different seasons in KMC areas. The main objectives of the study were: (1) To analyze the spatial variations of LST in the different wards as well as different LULC types during the different seasons; (2) to explore the optimal multiple regression model for each category of variables and the integration of all categories of variables; and (3) to identify the dominant influencing factor of each category of variables and among all categories of variables during the different seasons.

2. Materials and Methods

2.1. Descriptions of Study Area

Kolkata Municipal Corporation (KMC) is located on the east bank of the Hoogly River (Figure 1) in the state of West Bengal, India. It extends 22°27′12″ N to 22°37′56″ N and 88°14′33″ E to 88°27′32″ E, and covers an area of approximately 185.52 sq. km. The KMC contains 141 wards, and these small wards were chosen as a spatial unit for the study because they are the small administrative units under the jurisdiction of KMC for urban planning and management.

It is the gateway city of the eastern part of India and is the capital city of West Bengal. After independence, a rapid rate of urban growth was observed in Kolkata City and the surrounding areas due to the adoption of a mixed economy [36]. Originating as a colonial city, Kolkata City had a population of 2.92 million in 1961. However, according to the census of India 2011, it is the third most populous metropolitan city (4.49 million) in India with a population density of 24,252 per sq. km, and the population is expected to reach 4.97 million by 2022. (https://www.indiacensus.net/district/kolkata (accessed on 23 August 2022)). The rapid rate of urbanization has led to the transformation of natural land into artificial land [33], and as a result, the land surface temperature is increasing in the city and its surrounding areas [34]. As reported by [34], the land surface temperature has increased by 4.72 °C between 1990 and 2020 as the built-up area increased by 33.758% and the vegetation cover decreased by 25.521% in Kolkata and its adjacent areas.

According to the Köppen–Geiger climate classification, the climate is tropical wet and dry climate [37]. The average annual rainfall is 1655 mm and ranges between 11 mm (December) to 365 mm (July). The average annual temperature is 25.99 °C with a range of 19.1 °C (January) to 30.6 °C (May). Topographically, the city is near sea level with an average elevation of 6 m from the mean sea level. KMC has high elevation in the middle-northern parts of the study area, and approximately 98% of the study area has a height of 1 to 16 m. The elevation steadily decreases from the middle and middle-north towards the southwest, northeast, east, and southeast directions (Figure S1). The lowest elevation is observed in the eastern parts with marshes and swamplands.

2.2. Delineation of Seasons

KMC is located in a tropical wet and dry climate (AW) zone of India (Beck et al., 2018). There is seasonal variation in intensity and trend of LST in Kolkata and its adjacent areas [38]. Therefore, for the seasonal analysis of LST, KMC was delineated into three seasons, i.e., summer, transition, and winter. For delineating the seasons, monthly climate (temperature and precipitation) data from each year were obtained from meteorological stations (Behala airport station: At 1.80 m height from the ground) during the period of 1997–2019 (https://www.wunderground.com/ (accessed on 22 September 2021)). Thus, we collected 23 sets (23 years) of monthly temperature and precipitation data. Then, the average temperature and precipitation of each month were calculated based on the monthly data (23 sets of data) from 1997–2019. The average monthly temperature and precipitation were graphed to visualize as well as delineate the seasons. Figure 2 showed that the temperature in November, December, January, and February was low, with a temperature of <24 °C for these months. These four months also had the lowest precipitation—the precipitation was <30 mm. These months were the coldest and driest months of the year, thus these four months were considered as the winter season. Compared to the winter season, the temperature and precipitation in June, July, August, September, and October were higher, and the precipitation during these months was higher compared to other months. These months were the wettest months of the year, thus these five months were considered as the transition season. The temperature during April and May was higher compared to other months, and the precipitation was slightly higher than the winter season and lower than the transition season. Although the average temperature during March was lower than June, July, August, and September, the precipitation was much lower than these months (Figure 2) and the temperature was higher than the winter season. Therefore, the months of March, April, and May were the hottest months of the year, thus these three months were considered as the summer season.

2.3. Satellite Data and Descriptions

2.3.1. Landsat-8 Data

The USGS Landsat-8 surface reflectance product of 2019 was used for LULC classification and to estimate the biophysical indices, such as the Normalized Difference Build-up Index (NDBI), Enhanced Vegetation Index (EVI), and Modified Normalized Difference Water Index (MNDWI) in GEE platforms. The calibrated brightness temperature (band 10: Thermal infrared 1) of Landsat-8 was used to calculate LST (see Section 2.5. LST retrieval). These datasets were atmospherically corrected using LaSRC and include the cloud, shadow, water, and snow mask produced using CFMASK [39] from Landsat-8 OLI/TIRS sensors, which is available on the GEE repository. Landsat-8 has 11 spectral bands with 30 m spatial resolution and varying spectral resolutions (

Table 1. Spectral band properties of Landsat-8 OLI sensor.

Name	Description	Wavelength (nm)	Resolution (m)
B1	Coastal aerosol (CA)	433–453	30
B2	Blue (B)	450–515	30
B3	Green (G)	525–600	30
B4	Red (R)	630–680	30
B5	Near infrared (NIR)	845–885	30
B6	Sort wave infrared 1 (SWIR1)	1560–1660	30
B7	Sort wave infrared 2 (SWIR2)	2100–2300	30
B8	Panchromatic (PAN)	500–680	15
B9	Cirrus (C)	1360–1390	30
B10	Thermal infrared 1 (TIR1)	1030–1130	30
B11	Thermal infrared 2 (TIR2)	1150–1250	30

). Images with <10% cloud were considered, and we obtained a total of 77 images for this study.

2.3.2. LandScan Population Data

The population data were obtained from global population data released by the Oak Ridge National Laboratory (ORNL) of the US through the LandScan programme (https://landscan.ornl.gov/ (accessed on 22 September 2021)). This data provides an ambient population (average over 24 h) distribution at approximately 1 km (30″ × 30″) of spatial resolution [40]. LandScan was developed using the best available demographic (Census) and geographic data, remote sensing, and imagery analysis techniques within a multivariate dasymetric modelling framework to disaggregate the census counts within an administrative boundary [41]. We used the LandScan 2019 dataset to estimate the population in the study area.

2.3.3. VIIRS Nighttime Lights Data

The VIIRS nighttime version 1 data were used to estimate the nighttime light (NTL), which is available on GEE platforms; this is the monthly average radiance composite images that use nighttime data from the Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB). In 2011, NASA and NOAA launched the Suomi National Polar Partnership (SNPP) satellite, carrying the VIIRS, which collects low-light imaging data with 14-bit quantization, lower detection limits, a wider dynamic range, and approximately 45-times smaller pixels (463 m by 463 m) compared to the OLS [42]. In this study, the average DNB radiance (avg_rad) products of 2019 were used for the retrieval of the NTL in the study area.

2.3.4. OpenStreetMap Road Network Data

The road data were obtained from the open-source data website, OpenStreetMap (https://www.openstreetmap.org/ (accessed on 22 September 2021)). OpenStreetMap is a user-generated street map that contributes and maintains data about roads, trails, cafés, railway stations, and much more, all over the world [43]. In this study, the vector road data from OpenStreetMap in 2020 was merged to calculate the road density and was used for each season.

2.3.5. Sentinel-5P Air Pollutants Data

Sentinel-5P near real-time high-resolution imagery from the GEE repository was used to estimate the SO₂, NO₂, and CO. Sentinel-5P was launched on 13 October 2017 by the European Space Agency to monitor air pollution. The sensor carries a Tropospheric Monitoring Instrument (TROPOMI), which collects data daily with global coverage. Compared with the other air-quality observation satellite instruments, TROPOMI has a relatively high spatial resolution (1113 m × 1113 m) that is necessary for air quality applications [44]. Since the launch of Sentinel-5P, the TROPOMI dataset has been widely used and validated using in situ measurements with a low negative bias [45,46,47]. The Sentinel-5P products from 2019 were used for measuring the atmospheric pollutants.

The TROPOMI NO₂ processing system has been adopted from the algorithm developments for the DOMINO-2 product [48] and for the QA4ECV reprocessed dataset [49]. TROPOMI NO₂ represents the concentrations of collective nitrogen oxides, since during the daytime, i.e., in the presence of sunlight, a photochemical cycle involving ozone (O₃) converts NO into NO₂ and vice versa on a timescale of minutes (http://www.tropomi.eu/data-products (accessed on 17 October 2021)). In this study, the total vertical NO₂ column density product was used to estimate the NO₂ concentrations.

TROPOMI on the Sentinel-5 Precursor (S5P) satellite observes the CO global abundance by exploiting clear-sky and cloudy-sky earth radiance measurements in the 2.3 µm spectral range of the shortwave infrared (SWIR) part of the solar spectrum (http://www.tropomi.eu/data-products (accessed on 17 October 2021)). The operational processing for retrieving the total column density of carbon monoxide (CO) while simultaneously interfering with trace gases and effective cloud parameters (cloud height and optical thickness) was performed by the shortwave infrared carbon monoxide retrieval (SICOR) algorithm [50]. For retrieving the CO concentrations, the vertically-integrated CO column density product was used.

The TROPOMI SO₂ retrieval algorithm is based on the differential optical absorption spectroscopy (DOAS) technique [51]. For this paper, the SO₂ vertical column density at ground level was used to retrieve the SO₂ concentrations.

2.3.6. MCD19A2 Aerosol Optical Depth (AOD) Data

The MCD19A2 product was used to estimate the aerosol optical depth. The MCD19A2 is a version 6, level 2 Terra and Aqua combined Multi-Angle Implementation of Atmospheric Correction (MAIAC) Land Aerosol Optical Depth (AOD) product at 1 km of spatial resolution from MODIS daily observations. MAIAC is a new, advanced algorithm that uses time-series analysis and a combination of pixel- and image-based processing to improve the accuracy of cloud detection, aerosol retrievals, and atmospheric correction [52], and provide higher accuracy than MODIS dark target and blue algorithm over dark surfaces and smoke plumes [53]. The MCD19A2 products of 2019 were obtained from the GEE repository, and the green band of 0.55 μm of aerosol optical depth (AOD) product of MCD19A2 was used in this study because of its superior accuracy, as shown by previous studies [52].

2.3.7. SRTM Digital Elevation Data

Shuttle Radar Topography Mission (SRTM) digital elevation data that are available in GEE as USGS/SRTMGL1_003 were used to obtain topographical factors, such as slope, elevation, aspect, and hill shade. This SRTM DEM is a version 3 product at a resolution of 1 arc-second, which includes the topographic data from non-SRTM sources to fill the gaps (“voids”) in earlier versions of SRTM data [54]. SRTM digital elevation data are only available in 2000, therefore, SRTM digital elevation data from the 2000 launched mission were used.

2.3.8. MCD43A3 Albedo Data

The MCD43A3 data from 2019 was acquired from the GEE repository to measure the albedo. The MCD43A3 is a version 6, level 3 dataset from a daily 16-day product at 500 m of spatial resolution [55]. The MCD43A3 product combines registered, multidate, multiband, and atmospherically-corrected surface reflectance data from the MODIS and MISR instruments to fit a Bidirectional Reflectance Distribution Function (BRDF) in each spectral band [56]. The MCD43A3 Albedo Model dataset provides both directional hemispherical reflectance (black sky albedo) and bihemispherical reflectance (white sky albedo) for each of the MODIS surface reflectance bands (band 1 to band 7), as well as 3 broad-spectrum bands (visible, near-infrared, and shortwave). While the total energy reflected by the Earth’s surface in the shortwave domain is characterized by the shortwave (0.3–5.0 µm) broadband albedo, the visible (0.3–0.7 µm) and near-infrared (0.7–5.0 µm) broadband albedos are often also of interest due to the marked difference of the reflectance of vegetation in these two spectral regions. As black sky albedo (BSA) and white sky albedo (WSA) products of MCD43A3 are closely related [57], white sky albedo for shortwave was used in this study.

2.4. Land Use Land Cover Classification

The surface reflectance product from 2019 was used for classification, and the classification was conducted using cloud computing API in GEE platforms [58]. The annual mean composition was applied for each band because multi-image composition leads to more stable and accurate predictions than single date imagery [59]. While using individual bands (B2–B7), NDVI [60], NDBI [61], and MNDWI [62] were also performed for LULC classification. The topographic data from DEM was also included to distinguish between bodies of water and other types of LULC [63].

In this study, we used the random forest algorithm [64] for LULC classification in the GEE platform. The RF model has two parameters that are important for tuning, such as ntree and mtry [65]. Mtry is the number of predictor variables to be tested at each node, and the default value of the square root of the total predictors provides better results [64]. The second tuning parameter is ntree, which controls the total number of independent trees. Lots of trees are conducive to stabilizing the variable importance and variable interaction [65]. Therefore, we used the default value of mtry and 1000 ntree for the classification.

A total of 661 sampling points were visually delineated from very high-resolution Google Earth imagery (resolution < 1 m) using a stratified random sampling approach. These sampling points were collected based on the judgment of each LULC type, with consideration of no less than 50 sampling points per LULC category. Then, each sampling point was buffered by 30 m by 30 m to increase the reference data, as well as to enhance the classification result. After that, 70% of the data were used as training data to improve the classification, and the remaining 30% of the data were used as a validation dataset for accuracy assessment. Finally, we classified the study area into six types: artificial land, vegetation land, grassland, cropland, bodies of water, and barren land (

Table 2. Description of LULC types used in the classification.

LULC Types	Description
Artificial land (AL)	Commercial, residential, transport, and other socio-economic built-up areas
Vegetation land (VL)	Residential, recreational, aquatic, and roadside vegetation including shrubs and trees
Grass land (GL)	Open space with grass cover, playground, and other organized or unorganized grass cover areas
Crop land (CL)	Agricultural lands and fallow lands
Bodies of water (WB)	Surface covered by water including rivers, streams, ponds, lakes, canals, reservoirs, check dams, etc., which are considered as bodies of water
Bare land (BL)	All other lands including sandy or stony areas, dumping grounds, exposed soil surface, etc.

). The user’s accuracy (UA), producer’s accuracy (PA), F1 score, overall accuracy (OA), and Kappa’s coefficient (K) were calculated to assess the accuracy [66].

$$Kappa’s coefficient \left(K\right)=\frac{Overall accuracy-Agreement chance}{1-Agreement chance}$$

(1)

$$F1 score=2\times \frac{\mathrm{User}’\mathrm{s} \mathrm{accuracy}\times {\mathrm{Producer}}^{’}\mathrm{s} \mathrm{accuracy}}{\mathrm{User}’\mathrm{s} \mathrm{accuracy}+{\mathrm{Producer}}^{’}\mathrm{s} \mathrm{accuracy}}$$

(2)

2.5. LST Retrieval

The seasonal mean composition was applied for each band of Landsat-8 images to determine the LST for each season. Mean composition was applied because multi-image composition leads to more stable and accurate predictions than single date imagery [59]. The entire LST calculation was performed in the Google Earth Engine (GEE) platform. As GEE directly provides the brightness temperature (band 10) of Landsat-8, the LST was obtained across different seasons using Equation (1) [67]:

$$\mathrm{LST}=\frac{{\mathrm{T}}_{\mathrm{B}}}{1+\left(\frac{{\mathsf{\lambda } \mathrm{T}}_{\mathrm{B}}}{\mathsf{\rho }}\right)\ln \mathsf{\epsilon } }$$

(3)

where LST is the land surface temperature (LST) in Kelvin; λ is the wavelength of emitted radiance in meters (for which the peak response and the average of the limiting wavelengths $\mathsf{\lambda }$ = 11.5 μm) is used;${ \mathrm{T}}_{\mathrm{B}}$ is the brightness temperature in Kelvin (Landsat-8 band 10 brightness temperature was used); $\mathsf{\rho }$ = h.c/σ = 1.438 × 10⁻² mK, in which σ is the Boltzmann constant (1.38 × 10⁻²³ J/K); h is the Planck’s constant (6.626 × 10⁻³⁴ J s); and c is the velocity of light (2.998 × 10⁸ m/s) and ε is the emissivity.

Due to the complexity of the underlying surface of a city, it is difficult to determine the change of the emissivity value. Therefore, in this study, ε was simplified using the Equation (2) [68]:

ε = 1.0094 + 0.0047 ln NDVI

(4)

where NDVI was calculated from the red band and near-infrared band (Rouse et al., 1974). Furthermore, for pixels where NDVI < 0.157, emissivity was set as 0.923 [69]. For pixels where NDVI > 0.727, emissivity was set as 0.9925 [70]. The LST’s unit was converted to degrees Celsius using the relation of 0 °C equals 273.15 K. Then, the average LST °C value in each ward was calculated for each season.

2.6. LST Distributions

To analyze the ward-wise spatial distribution of LST, the LST data were divided into five categories (

Table 3. Five categories of LST values.

LST Category	Range of Each Level
Very low	$T \le T_{mean }-1.5std$
Low	$T_{mean }-1.5std<T \le T_{mean }-0.5std$
Moderate	$T_{mean }-0.5std<T \le T_{mean }+0.5std$
High	$T_{mean }+0.5std<T \le T_{mean }+1.5std$
Very high	$T>T_{mean }+1.5std$

) for each season based on the mean and standard deviation value [71]. For the LST distribution in different LULC types, the mean and standard deviation of LST of each LULC type was calculated for the different seasons. Furthermore, the average LST was calculated for the different seasons.

2.7. Influencing Factors Selection

A seasonal (summer, transition, and winter) mean composition was applied for Landsat-8, VIIRS nighttime light, Sentinel-5P air pollutants, MCD19A2 aerosol optical depth, and MCD43A3 albedo data to get seasonal information. As LandScan data are an annually published product, the LandScan 2019 dataset was used as the influencing factor for three seasons. SRTM digital elevation data from a launched mission from 2000 and OpenStreetMap road network data from 2020 were used as the influencing factors for three seasons. Then, a total of 25 indicators were determined for each season as influencing factors and divided into 6 categories, such as pollutant parameters, surface properties, topographical parameters, landscape compositions, landscape configurations, and socio-economic parameters (

Table 4. Influencing factors of land surface temperature (LST) used in this study.

Categories of Variables	Variables	Meaning of Variables
Pollutant parameters	AOD	Mean value of aerosol optical depth in the ward
	NO2	Mean value of NO2 concentration in the ward (mol/m2)
	CO	Mean value of CO concentration in the ward (mol/m2)
	SO2	Mean value of SO2 concentration in the ward (mol/m2)
Surface properties	Albedo	Mean value of white scale albedo in the ward
	EVI	Mean value of EVI in the ward
	MNDWI	Mean value of MNDWI in the ward
	NDBI	Mean value of NDBI in the ward
Topographical parameters	Slope	Mean value of slope in the ward (degree)
	Elevation	Mean value of elevation in the ward (meter)
	Hill shade	Mean value of hill shade in the ward (degree)
	Aspect	Mean value of aspect in the ward (degree)
Landscape configurations	PD	Density of the landscape patches in the ward (number/100 ha)
	ED	Edge density of the landscape patches in the ward (m/ha)
	LSI	Landscape shape index of the landscape in the ward
	MSI	Mean shape index of the landscape patches in the ward
	SHDI	Shannon diversity index of the landscape in the ward
Socio-economic parameter	RD	Road density in the ward (km/sq.km)
	NTL	Mean value of nighttime lights in the ward (nW/cm2/sr)
	Pop	Mean value of population in the ward
Landscape compositions	VL	Percentage of vegetation land in the ward
	GL	Percentage of grassland in the ward
	WB	Percentage of bodies of water in the ward
	BL	Percentage of bare land in the ward
	AL	Percentage of artificial land in the ward

). The average value of all the influencing factors for each ward were calculated for three seasons. Pixels with more than 50% of the area within the ward boundary were considered for calculation of the average value.

The urban surface cover properties are the most important influential factors of LST [72] because each surface cover has unique radiative, thermal, moisture, and aerodynamic properties [73]. EVI, MNDWI, NDBI, and albedo were selected as influencing factors of LST to cover the urban surface properties. Among surface properties, NDVI [60] is the most popular and widely-used index to understand the effect of vegetation cover on LST [22]. However, in this study, EVI [74] was used because the Normalized Difference Vegetation Index (NDVI) is prone to saturation [28]. The surface water cover information is the sensitive index for detecting LST change [75]. The most popular water detection index is NDWI [76], but the Normalized Difference Water Index (NDWI) is integrated with construction land features. To overcome the deficiency, MNDWI [62] was used to characterize the surface water. Previous studies showed that NDBI has a significant effect on LST [77]. Thus, NDBI was used to identify the relationship with LST across the cities. Albedo in the urban environment can also influence the LST [78], so white sky albedo (WSA) was used to explore the effect on LST. The seasonal average value of the biophysical indices of each ward were calculated for use as influencing factors of LST.

$$EVI=2.5\times \left(\frac{B5-B4}{B5+\left(6\times B4\right)-\left(7.5\times B2\right)+1}\right)$$

(5)

$$MNDWI=\frac{\left(B3-B5\right)}{\left(B3+B5\right)}$$

(6)

$$NDBI=\frac{\left(B6-B5\right)}{\left(B6+B5\right)}$$

(7)

Many previous studies have shown that different landscape compositions and landscape configuration indicators have a significant impact on urban LST [79,80]. Therefore, five commonly-used landscape compositions and five landscape configuration indicators were selected as influencing factors of LST. The five landscape composition variables are the percentage of vegetation land (VL), artificial land (AL), grassland (GL), bodies of water (WB), and bare land (BL). The five landscape configuration indicators are patch density (PD), edge density (ED), landscape shape index (LSI), mean shape index (MSI), and Shannan diversity index (SHDI). The landscape compositions and configurations were calculated for each ward for use as influencing factors of LST. The software, Fragstats 4.0, was used to calculate the landscape configurations from the LULC map.

A number of studies have also shown that different atmospheric pollutant parameters have had a significant impact on LST. Basically, the concentration aerosol optical depth (AOD), NO₂, SO₂, and CO can significantly influence the LST in built-up areas [30,31,32,81]. Therefore, in this study, AOD, NO₂, SO₂, and CO were selected and the seasonal average value of each parameter of each ward was calculated.

Different topographical parameters from different environments could significantly affect LST [25,29]. Therefore, four important topographical parameters, such as elevation, slope, aspect, and hill shade, were selected to analyze the influence of topographical factors on LST. After that, the average value of the topographical parameters of each ward were calculated to use as influencing factors of LST.

Different socio-economic factors such as nighttime lights (NTL), population (Pop), and road density (RD) have significant effects on LST [23,24]. Therefore, the present study also examined the effects of NTL, population, and road density on LST in urban environments. Then, each socio-economic factor was estimated for each ward for use as influencing factors.

2.8. Factors Analysis

To examine the effect of each factor on LST variation in each season, the most common and widely-used ordinary least square regression was used. For each season, the average LST of each ward was taken as a dependent variable, and 25 influencing factors were used as independent variables (Table 4). The significance test at the 5% level was used to identify significant influencing factors of LST using this method.

The multiple regression model was developed to model the LST based on the combination of explanatory variables. The major problem of the multiple regression model is noisy variables [82]. Therefore, we only selected variables that must have been statistically significant at a threshold of ≤5%. Although, it is possible to select optimal multiple regression models using the stepwise regression method, it is not able to consider all possible combinations with avoiding the collinearity of influencing factors [83]. To address this problem and take all possible combinations of influencing factors for more information, the “Leaps” package [84] was used in the R statistical software. This package helps regression subset selection using an efficient branch and bound algorithm. We used one best model for each number of explanatory variables (nbest = 1) with no limit on the number of variables (nvmax = NULL) using the “exhaustive” method. Then, an adjusted r-square was taken to identify the best subset of the model. Then, a multiple r-squared was calculated for the best model to determine the explanatory rate of the model. The “Leaps” package was applied individually for six categories, such as pollutant parameters, surface properties, topographical parameters, landscape compositions, landscape configurations, and socio-economic parameters, to explain the explanatory rate of each category of variables during different seasons. Finally, the “Leaps” package was applied for the summer, transition, and winter season using all influencing factors of LST.

After that, using the “hier.part” package in the R statistical software, hierarchical partitioning (HP) was conducted to identify the independent effects of explanatory variables on LST. HP is more suitable for environmental data analysis because the HP is not affected by the multicollinearity of the variables [85]. However, the HP function using “hier.part” produces a “minor rounding error” for models constructed with more than nine predictors [86], wherein datasets are likely to interchange their positions. To overcome this problem, we ran the model 100 times with the variables entered in a different order; then, the average of the 100 runs was taken to detect the independent effects.

3. Results and Discussion

3.1. LULC Analysis

The area and accuracy of each LULC type are presented in

Table 5. Area of each LULC type and the user’s, producer’s, F1 score, overall accuracy, and kappa coefficient for LULC classification.

	Area			Accuracy
LULC Types	Area (Km2)	Area (%)	User’s Accuracy (%)	Producer’s Accuracy (%)	F1 Score	Overall Accuracy	Kappa’s Coefficient (K)
Artificial land (AL)	148.66	80.13	94.82	96.83	95.82
Vegetation land (VL)	13.75	7.41	95.45	91.31	93.33
Grass land (GL)	15.16	8.17	84.61	83.02	83.81	92.64	0.87
Crop land (CL)	0.71	0.38	85.71	81.82	83.72
Bodies of water (WB)	6.32	3.41	93.33	90.32	91.80
Bare land (BL)	0.92	0.50	84	80.77	82.35

. We had good LULC classifications with an overall accuracy and Kappa’s coefficient of 92.64% and 0.87, respectively. Artificial land, vegetation land, and bodies of water have a higher producer’s accuracy, user’s accuracy and F1 score. Artificial land was the largest distributed land cover (80.13%), followed by grassland (8.17%), vegetation land (7.41%), bodies of water (3.41%), bare land (0.50%), and cropland (0.38%) (Figure 3 and Table 5).

3.2. Spatial Variation of LST in Relation to LULC Types

The rapid unplanned expansion of urbanization is a major issue in the urban areas of KMC [33], which is increasing the artificial area on a large scale, and as a result, the LST is increasing in Kolkata and its surrounding areas [34]. Therefore, information about the variation and pattern of LST in the different LULC types in every small spatial unit of a city is important on seasonal scales for local-scale adaptation and mitigation [87]. The highest LST corresponded to artificial land, followed by bare land, vegetation land, grassland, cropland, and bodies of water (Figure 4). The average LST distribution in the different LULC types followed similar patterns for the different seasons. The average LST of each LULC type decreased continuously from summer to transition, then to winter (Figure 4), which may be due to the decrease in solar radiation from summer to transition to winter because of the tilted axis of rotation of the Earth. The average LST was 30.432 °C in summer, followed by27.415 °C in transition, and 24.894 °C during the winter season. Artificial land corresponded to the highest LST in all seasons because built-up lands have different thermal bulk properties (heat capacity and thermal conductivity) and radiative properties (albedo and emissivity) [88]. For this reason, the urban artificial surface absorbs significantly more solar radiation, which causes a change in the energy budget of the urban area, often leading to a higher LST [89]. Dense buildings also reduce wind speed and block the longwave emissions from the ground, which helps to gather heat inside the city and increase the temperature [90]. The second-highest LST was observed in bare lands. This result is quite similar to the study in [90], wherein they revealed that bare soil’s heat capacity is minimal with a lower moisture content and evapotranspiration rate, resulting in the rapid rate of heating during the daytime. Therefore, high temperature surfaces were developed in exposed areas. The lowest temperatures were observed in bodies of water, followed by croplands, grasslands, and vegetation lands. Plants and bodies of water act as an evaporation cooling system by spreading the moisture to the atmosphere through evaporation, resulting in a reduction of the LST [90]. Furthermore, due to the higher heat capacity of bodies of water than land surface, water surface temperatures increase slowly, so low temperatures are observed in bodies of water compared to land surfaces under the same solar radiation [90]. Vegetation also creates an albedo effect that is 15% higher than urban surfaces due to the small absorption of heat and higher reflectance [91]. Furthermore, the canopy of vegetation increases the shade effect, which reduces incident radiation [89]. The average temperature of croplands was lower than that of grasslands and vegetation lands, which may be due to the croplands acting as vegetation and irrigation lands, and irrigation decreases land surface temperatures during daytime [92]. Therefore, green surfaces and bodies of water are the major factors that mitigate the increase in LST. This result is quite similar to the study from [93], which revealed that high temperatures and low temperatures correspond to built-up areas and bodies of water in three megacities in China.

Therefore, it is clear from Figure 5 and Figure 6 that the high-temperature ward as well as the regions were mostly distributed in the middle-northern parts and somewhere in the western parts of highly urbanized areas, where artificial land cover predominated, lacked bodies of water and green cover (Figure 3). There is a large number of previous studies [77,93] that also revealed that artificial land surfaces are highly related to a high LST. The number of very high temperature wards was 5, 12, and 6 in the summer, transition, and winter seasons, respectively (Figure 6). The number of very low temperature wards was 9, 3, and 11 in the summer, transition, and winter seasons, respectively (Figure 6). Low and very low temperature wards were distributed in the surrounding areas that dominate the eastern part of the study area, which is most likely due to the distribution of green zones and bodies of water (Figure 3).

3.3. Effects of Each Category of Influencing Factors in the Different Season

The OLS regression analysis showed that many influencing factors were significantly fitted across the cities for the different seasons (

Table 6. R-squared of ordinary least square regression (OLS) for influencing factors affecting LST.

	Influencing Factors	Summer	Transition	Winter
Pollutant parameters	AOD	0.368 ***(−)	0.444 ***(−)	0.315 ***(−)
	NO2	0.635 ***(+)	0.423 ***(+)	0.004(−)
	CO	0.575 ***(+)	0.005(+)	0.055 *(+)
	SO2	0.398 ***(−)	0.208 ***(−)	0.165 ***(−)
Surface properties	Albedo	0.202 ***(−)	0.316 ***(−)	0.341 ***(−)
	EVI	0.689 ***(−)	0.569 ***(−)	0.708 ***(−)
	MNDWI	0.059 **(−)	0.025(−)	0.072 **(−)
	NDBI	0.765 ***(+)	0.671 ***(+)	0.653 ***(+)
Topographical parameters	Slope	0.081 ***(+)	0.149 ***(+)	0.138 ***(+)
	Elevation	0.294 ***(+)	0.299 ***(+)	0.228 ***(+)
	Hill shade	0.013(−)	0.016(−)	0.015(−)
	Aspect	0.034 *(+)	0.031*(+)	0.030 *(+)
Landscape configurations	PD	0.621 ***(−)	0.551 ***(−)	0.494 ***(−)
	ED	0.556 ***(−)	0.520 ***(−)	0.448 ***(−)
	LSI	0.559 ***(−)	0.524 ***(−)	0.452 ***(−)
	MSI	0.000(−)	0.007(−)	0.001(−)
	SHDI	0.587 ***(−)	0.521 ***(−)	0.449 ***(−)
Socio-economic parameters	RD	0.210 ***(+)	0.202 ***(+)	0.221 ***(+)
	NTL	0.385 ***(+)	0.411 ***(+)	0.441 ***(+)
	Pop	0.572 ***(+)	0.531 ***(+)	0.636 ***(+)
Landscape compositions	AL	0.571 ***(+)	0.397 ***(+)	0.477 ***(+)
	VL	0.374 ***(−)	0.261 ***(−)	0.274 ***(−)
	GL	0.438 ***(−)	0.310 ***(−)	0.437 ***(−)
	WB	0.191 ***(−)	0.117 ***(−)	0.119 ***(−)
	BL	0.022(−)	0.002(−)	0.019(−)

***, **, and * represent the significance at 0.001, 0.01, and 0.05 levels, respectively. (+) and (−) indicates that the correlation between the influencing factors and LST is positive and negative, respectively. Bold color indicates the most significant factors of each category of influencing factors.

). The relationship between the LST and influencing factors during the different seasons are shown in Figure S2. The significantly fitted or correlated variables are those that are statistically significant at the 5% significance level and are later chosen for all-subsets regression modelling.

The surface properties characterize the highest explanatory rate compared to the other categories of variables for all seasons (

Table 7. Results of optimum regression model for each category of variables in different seasons.

Summer Season
Categories of Variables	Model	Adj. R2	Mul. R2
Surface properties	Y = 33.9511 − 0.8384 × EVI + 12.0250 × NDBI − 13.2100 × Albedo	0.785	0.789
Pollutant parameters	Y = −10.9900 + 61324.1600 × NO2 + 694.6000 × CO	0.644	0.649
Topographic parameters	Y = 28.9178 + 0.2326 × Elevation	0.289	0.294
Landscape configurations	Y = 34.1219−0.0339 × PD + 0.0087 × ED	0.641	0.645
Landscape compositions	Y = 26.4537 + 0.0525 × AL	0.568	0.571
Socio-economic parameters	Y = 29.0500 + 0.00002 × Pop + 0.0310 × NTL	0.625	0.632
Transition season
Categories of variables	Model	Adj. R2	Mul. R2
Surface properties	Y = 33.1890 + 14.2030 × NDBI − 36.5320 × Albedo	0.701	0.705
Pollutant parameters	Y = 18.0900 + 01 + 111500.0127 × NO2 − 8172.1540 × SO2 − 2.7741 × AOD	0.597	0.606
Topographic parameters	Y = 25.0790 + 0.2394 × Elevation + 0.3174 × Slope	0.304	0.314
Landscape configurations	Y = 31.3880 − 0.0359 × PD + 3.2742 × SHDI	0.554	0.561
Landscape compositions	Y = 23.5827 + 0.0515 × AL	0.393	0.398
Socio-economic parameters	Y = 26.2710 + 0.00002 × Pop + 0.0177 × NTL	0.548	0.556
Winter season
Categories of variables	Model	Adj. R2	Mul. R2
Surface properties	Y = 30.2312 − 01.6033 × EVI + 6.2838 × NDBI − 27.9831 × Albedo + 0.1589 × MNDWI	0.779	0.786
Pollutant parameters	Y = 37.6200 − 0.0130 × AOD − 8703.0312 × SO2	0.503	0.511
Topographic parameters	Y = 23.1459 + 0.1635 × Elevation + 0.3187 × Slope	0.239	0.251
Landscape configurations	Y = 42.5701 − 0.0464 × PD + 0.0282 × ED + 5.2900 × SHDI − 10.8610 × LSI	0.53	0.543
Landscape compositions	Y = 26.0201 − 0.0629 × GL − 0.0353 × VL − 0.0299 × WB	0.483	0.495
Socio-economic parameters	Y = 23.5302 + 0.00002 × Population + 0.0197 × NTL	0.677	0.682

). Various studies also used surface biophysical properties to explain the variability of LST [24,26], and reported that the surface cover properties are the most influencing factors of urban LST [72] because the surface biophysical indices contain more information about particular properties [24]. Our results are also consistent with these studies. The explanatory rates (R-squared) of the surface properties were 78.9%, 78.6%, and 70.5% in the summer, winter, and transition seasons, respectively. This indicates, along with increasing rainfall, that the explanatory rate of surface properties has decreased, resulting in surface properties being more important for explaining the LST in the winter and summer season when the temperature is relatively low and high with low rainfall. This result agreed quite well with the study from [24], which reported that biophysical parameters could explain more of the variation in LST during low- or high-temperature seasons. In the summer season, the combination of NDBI, EVI, and albedo gave the highest explanatory rate (Figure 7A), and NDBI was the most significant contributing factor (50%), followed by EVI (40%) and albedo (10%). In the transition season, the combination of NDBI and albedo gave the highest explanatory rates with an independent effect of 73% and 27%, respectively (Figure 7B). In the winter season, the highest explanatory rate was received from the combination of all surface properties, and NDBI was the most significant contributing factor (41%) compared to EVI (39%), albedo (15%), and MNDWI (5%) (Figure 7C). Compared to the summer and transition seasons, the number of factors was increased to explain the variation in LST during the winter season, which indicates that during low temperatures, more factors were required to explain the variation in LST. Conversely, in the transition season, only a few variables of surface properties (NDBI and albedo) are entered into the optimal model to explain the variation in LST with an explanatory rate of 70.5%, which indicates that there are other robust factors that control the variability of LST. NDBI and albedo were constantly selected for the optimal model for all the seasons, which indicates that the combination of NDBI and albedo is important to explain the more variation of LST for all seasons. EVI was selected for the summer and winter seasons, and MNDWI was selected for the winter—in combination with other variables—which indicates that EVI and MNDWI have seasonal importance for explaining the LST in combination with other surface properties. Furthermore, the OLS regression explanatory rates of EVI and NDBI were higher with a high significance (Table 6), which is consistent with previous studies. Various studies also reported that MNDWI or bodies of water have a highly significant effect on LST [26,90], but in this study, MNDWI had a very low explanatory rate, which was mainly due to the total area of bodies of water being small (3.41%) and concentrated in one particular site (eastern parts) of the study area (Table 5 and Figure 3).

In the categories of pollutant parameters, different factors were entered for the optimal models for the different seasons. Previous studies have always ignored the category of pollutant parameters, but pollutant parameters could explain the 64.9% variation in the LST during the summer, 60.6% variation during the transition, and 51.1% variation during the winter season (Table 7). This indicates that the explanatory rate of pollutant parameters has decreased with the decrease in LST, thereby resulting pollutant parameters could explain more variation of LST when the temperature is relatively high in the study area. In the summer season, only CO and NO₂ were entered into the optimal model with an independent effect of 55% and 45%, respectively (Figure 7A). In the transition season, the combination of NO₂, AOD, and SO₂ gave the highest explanatory rate, and NO₂ had the largest independent effect (44%), followed by AOD (37%) and SO₂ (19%) (Figure 7B). In the winter season, AOD and SO₂ gave the highest explanatory rate, with an independent effect of 65% and 35%, respectively (Figure 7C). Therefore, more factors were required to explain the variation in LST during the rainy season, and limited factors were required during the relatively low or high temperatures of low rainfall seasons. That is to say, AOD and SO₂ could explain the high levels of variation in LST when the temperature was relatively low in the study area. AOD and SO₂ were also the highest explanatory factors during the winter season in OLS regression (Table 6). However, when the summer season temperatures are high, the combination of NO₂ and CO could explain the high levels of variation in LST. Furthermore, in the OLS regression, NO₂ was the highest explanatory factor, followed by CO, SO₂, and AOD (Table 6) during the summer season.

Regarding socio-economic factors, the explanatory rates were greater in winter (68.2%), followed by the summer (64.9%) and transition season (60.6%) (Table 7). The combination of population and NTL gave the largest explanatory rates for all seasons with independent effects of 66% and 34% in the summer (Figure 7A), 62% and 38% in the transition (Figure 7B), and 68% and 32% in the winter (Figure 7C), respectively. Furthermore, in OLS regression, population and NTL were strongly correlated with LST for all seasons (Table 6), which is consistent with previous studies [23]. The explanatory rate of road density was 20.99%, 20.25%, and 22.14% during the summer, transition, and winter seasons, respectively (Table 6). The higher explanatory rate in winter indicates that when LST is relatively low (winter season), the effect of socio-economic parameters on LST variation is high. Various prior studies also reported that the amount of anthropogenic heat released over cities due to socio-economic activities is greater in winter than in summer [94,95], which indicated that regulating the population size could mitigate the LST [96]. Basically, nighttime lights are always associated with a high per capita income [97], which indicates high economic activity [98].

The explanatory rate of LST variation by topographic factors was relatively low compared to other groups of variables (Table 7). Only the elevation factor entered the optimal model for the summer season, which was also the higher interpreter compared to other topographic factors in the OLS regression (Figure 7A). For the transition and winter seasons, combinations of elevation and slope gave the highest explanatory rate with an independent effect of 74% and 26% in the transition, and 68% and 32% in the winter, respectively (Figure 7B,C); the explanatory rate was higher in transition (31.4%) than in summer (29.4%) and winter (25.1%). However, all the significant topographic factors were positively related with LST (Table 6), which is opposite to the study by Peng et al. (2020) that reported that topographic parameters, such as slope and elevation, generally have a negative effect on LST because as elevation increases, the air becomes thinner and loses its capability to retain heat. However, our study depicts opposite results, and this was mainly because of two reasons: firstly, a very small variation in topography (2.42 to 17.48 m) throughout the study regions, so there may be no significant difference in air density in the study region; secondly, the low elevation and low slope areas were covered by green covers and bodies of water, and the high elevation and high slope areas were covered by human activities as well as artificial lands. Therefore, artificial lands with a high elevation and sloping regions had a negligible infiltration rate and helped to quickly dry up quickly downward water flows due to the force of gravity. Given this, the lower availability of moisture helped to quickly heat up elevated and sloping regions, and the high availability of moisture helped to reduce the temperature in low land areas. This might also cause the topographic parameters that are largely explained by the variability in LST during the transition season compared to the summer and winter seasons.

Landscape configurations have higher explanatory power compared to landscape compositions, which explain the variation in LST for all seasons (Table 7). Several previous studies [25,99] have reported that landscape configurations could explain more of the variation in LST than landscape composition. Again, [24,26] reported that the impact of landscape configurations is limited compared to landscape compositions. Therefore, it is not possible to make any reasonable comparisons with other studies in different cities as different regions and study periods may differentiate the effect of landscape configurations on LST. Except for MSI, all other configuration factors had a strong negative correlation with LST in OLS regression (Table 6). The explanatory rates of the optimum model were 64.5%, 56.1%, and 54.3% for the summer, transition, and winter, respectively. This also indicates, along with the decrease in LST, that the explanatory rate of the landscape configurations decreased, and as a result, the landscape configurations could better explain the higher level of variation in LST when the temperature was high in the study area. During the summer season, the combination of PD, ED, and SHDI gave the highest explanatory rate (Figure 7A), and the PD was the largest contributor (38%), followed by SHDI (32%) and ED (30%). During the transition season, we found that the highest explanatory rate was from the combination of PD and SHDI, which have an independent effect of 53% and 47%, respectively (Figure 7B). During the winter season, the combination of all significant landscape configurations, such as PD, ED, SHDI, and LSI, gave the highest explanatory rates, with independent effects of 33%, 21%, 24%, and 22%, respectively (Figure 7C). It is also noted that in the winter season, the highest number of factors were required of landscape configurations for the optimal model with the lowest explanatory rate compared to other seasons to explain the variation in LST, thereby indicating that the variation in LST became more complicated during the winter season. PD and SHDI were entered into the optimal model for all seasons, given that the combination of PD and SHDI was important for explaining the high levels of variation in LST for all seasons. However, ED was entered into the optimal model for the summer and winter seasons, and LSI was only entered for the winter, which indicates that ED and LSI have seasonal importance that explains the higher levels of variation in LST in combination with the other landscape configurations.

The explanatory rates for landscape compositions were 57.1%, 39.8%, and 49.5% during the summer, transition, and winter, respectively (Table 7). In the summer and transition season, only artificial land (AL) entered into the optimal model (Figure 7A,B), but in winter, the combination of vegetation land (VL), grassland (GL) and bodies of water (WB) gave the highest explanatory rates with an independent effect of 29%, 61%, and 10%, respectively (Figure 7C). This result indicates that the variation in LST was more likely affected by one dominant landscape composition factor (artificial land) when the temperature was relatively high, but when the temperature was low, the proportion of green space cover and bodies of water could explain more of the variation in LST, thereby resulting in the higher number of variables were required to explain the variation in LST during the winter season. The artificial land in OLS regression also gave the highest explanation in summer than in other seasons (Table 6). This finding is consistent with previous studies, which suggest that the higher the proportion of artificial land, the stronger the effect on LST during the summer [24] due to the intense heating of urban artificial surfaces [100,101]. During the winter season, there is a lack of adequate solar radiation, so the combination of green cover areas and bodies of water play a greater role compared to the single-level artificial surface that has been used to explain the variation in LST. Although bodies of water and bare lands were associated with a lower LST and higher LST, respectively (Figure 4), the explanatory rates of bodies of water were only 19.09%, 11.75%, and 11.91% during the summer, transition, and winter, respectively. Moreover, the explanatory rate of bare land was not significant for any of the seasons. This may be due to the low availability of bodies of water and bare lands in the study area (Figure 3 and Table 5).

However, in the summer and transition season, surface properties gave highest explanatory rates, followed by pollutant parameters, landscape configurations, socio-economic parameters, and topographic parameters. In contrast, in the winter season, surface properties gave the highest explanatory rates, followed by socio-economic parameters, pollutant parameters, landscape configurations, and topographic parameters. In addition, the explanatory rate (R-squared value) of all factors in the OLS regression varied for the different seasons. Therefore, the effects of influencing factors on LST are seasonally changed.

3.4. Integrated Effects of Al l Influencing Factors in the Different Seasons

In heterogeneous urban environments, LST is not controlled by a single factor or by single categories of factors [102]. Previous studies [24,26] used a large number of influencing factors to explain the variability in LST, but in this study, we have additionally used two more categories of influencing factors, pollutant parameters and topographic parameters. A total of 25 influencing factors from 6 categories of variables were used in this study, and 22, 20, and 21 influencing factors were included for all-subsets regression for the summer, transition, and winter season, respectively.

Ten significant factors, NO₂, albedo, EVI, MNDWI, NDBI, slope, LSI, NTL, population, and VL, were entered into the optimal model for the summer season (Figure 8A), for which the explanatory rate was 89.4% (

Table 8. Results from optimum regression model of integration of all influencing factors for different seasons.

Summer Season
Model Coefficient	Value	Adjusted R-Squared	Multiple R-Squared
Intercept	27.718	0.885	0.894
NO2	3515.894
Albedo	−1.920
NDBI	6.641
VL	−0.001
MNDWI	−0.694
EVI	−0.578
Slope	−0.103
LSI	−0.405
NTL	0.016
Pop	0.000
Transition Season
Model coefficient	Value	Adjusted R-squared	Multiple R-squared
Intercept	45.153	0.795	0.814
SO2	−4325.929
NDBI	17.909
Elevation	−0.061
Slope	0.533
PD	0.007
ED	0.026
LSI	−12.203
NTL	0.007
Pop	0.000
AL	−0.019
GL	−0.010
Winter Season
Model coefficient	Value	Adjusted R-squared	Multiple R-squared
Intercept	31.250	0.874	0.887
AOD	−1.897
SO2	−2314.895
NDBI	9.718
MNDWI	−1.270
EVI	−1.171
Elevation	−0.052
LSI	−2.493
SHDI	2.798
NTL	0.013
Pop	0.000
AL	−0.006
WB	0.014

). For the transition season, 11 explanatory variables were included, SO₂, NDBI, slope, elevation, PD, ED, LSI, NTL, population, AL, and GL (Figure 8B), for which the explanatory rate was 81.4% (Table 8). For the winter season, the explanatory rate was 88.7% (Table 8) and 12 variables were included: AOD, SO₂, EVI, MNDWI, NDBI, elevation, LSI, SHDI, NTL, population, AL, and WB (Figure 8C). Therefore, from summer to transition to winter, the number of factors for the optimal model was increased, indicating that with a decrease in the LST, more factors were required to explain the variation in LST. This result is consistent with the study from [24], which revealed that more factors are required to explain the variability of LST from the warm to the cold season. For the transition season, the explanatory rate was far lower (81.4%) than the summer (89.4%) and winter (88.7%), which indicates that there are other important factors controlling the variation in LST during rainy periods. As more than 9 influencing factors were entered into the optimal models for all seasons, we ran the HP model 100 times with different orders of the set of variables, as recommended by [86], then the average of the 100 runs was taken to identify the independent effects. The independent effects of the optimum models for each season are shown in Figure 9. NDBI was the dominant contributing factor (19.49%), followed by EVI (16.23%), NO₂ (15.22%), population (13.15%), LSI (11.79%), VL (9.78%), NTL (8.12%), albedo (3.69%), slope (1.59%), and MNDWI (0.93%) for the summer. Therefore, in the summer season, the sum of surface properties independently contributed 39.41% (NDBI, EVI, and albedo), followed by socio-economic parameters (21.27%), pollutant parameters (15.22%), landscape configuration (11.79%) and compositions (10.71%), and topographic parameters (1.59%). For the transition season, NDBI also was the dominant factor (15.95%) for LST variations, while PD, population, LSI, ED, NTL, AL, GL, elevation, SO₂, and slope contributed 11.16%, 10.81%, 10.53%, 10.31%, 8.58%, 8.27%, 7.33%, 6.29%, 5.87% and 4.93%, respectively. For the transition season, the sum of landscape configuration made the largest contribution (32%), followed by socio-economic parameters (19.39%), surface properties (15.95%), landscape compositions (15.6%), topographic parameters (11.22%), and pollutant parameters (5.87%). For the winter season, EVI and NDBI were the dominant contributors (14.78% and 14.66%, respectively), whereas the independent effects of the other selected factors were 12.31%, 9.26%, 9.13%, 8.47%, 8.37%, 6.42%, 4.74%, 4.37%, 4.27%, and 3.22% for population, AL, NTL, LSI, SHDI, AOD, elevation, SO₂, WB, and MNDWI, respectively. In the winter season, the sum of surface properties contributed 33.66%, followed by socio-economic parameters (21.44%), landscape configurations (16.84%) and compositions (13.63%), pollutant parameters (10.79%), and topographic parameters (4.74%). Therefore, surface properties largely contributed during the summer season, compared to the winter season when the temperature was relatively high and low with less rainfall. However, landscape configurations largely contributed during the transition season when the temperature was medium with high rainfall. The socio-economic properties (19.39–21.44%) and landscape compositions (10.71–15.6%) had slightly consistent contributions for all the seasons. The pollutant parameters and topographic parameters also largely contributed to LST during the summer (15.22%) and transition seasons (11.22%), respectively. Therefore, there were significant differences for the contributions from the influencing factors for different seasons. This result is similar to the study from [24,103] and revealed that various influencing factors behave differently for the different seasons. The first five factors contributed 75.89%, 58.76%, and 60.14% for the summer, transition and winter seasons, respectively. This result also indicated the high temperature (summer season) variation in LST was more likely affected by the combination of a few dominant factors, and during the transition and winter season, the independent effect tended to spread among the variables.

4. Limitations and Future of LST Studies

Although we have used a number of influencing factors to explain the variability in LST for different seasons, there are still a few limitations. The inclusion of other factors, such as climatic conditions [104], urban morphology [105], energy consumption rate [106], and distance to heat source and sink [107], could make future work more robust and informative. In this study, small administrative units were taken as spatial units. However, the size of every spatial unit was not equal, and each unit had unique features that could affect LST [108]. Subsequent research should focus on LST under different grid sizes (resolution) with a large number of influencing factors. Further research is also needed on different statistical methods for factor analysis. Different non-linear and non-parametric models should be used for their higher flexibility and performance compared to linear models, and comparisons between the models could help to discover the best model. It would also be interesting to conduct similar analyses in different spatial–temporal contexts to gain insights into evolution.

5. Implications for Managements

Changes in LST have become an important issue in urban ecosystems [109] as well as on human beings [110]. We can’t stop the process of urbanization, so we should adopt appropriate adaptation and mitigation strategies to balance the environment and LST. To adopt adaptation and mitigation strategies, a few general pieces of information are required: (1) What are the variations and patterns in LST for different seasons; and (2) what are the significant controlling factors of urban LST for different seasons? This study provides detailed information about the spatial variation and pattern of LST related to different LULC types for each small administrative unit of KMC during different seasons. This study also depicts a better understanding of the effects of categories of influencing factors on LST during the different seasons. This information is important for undertaking adaptation and mitigation strategies in urban areas via providing guidelines for governments, NGOs, and stakeholders. Various studies show the spatial patterns of LST in relation to a single influencing factor or few limited categories of influencing factors. To overcome the previous limitations, this study used 25 influencing factors from 6 categories of variables and used OLS and all-subsets regression to find out the significant factors and optimum model of LST. In addition, to avoid the “minor rounding error” for determining the independent effect of each factor, we ran the model 100 times with different orders. Thus, we have determined the effects of each category of variables and integrated effects of all variables on LST for the different seasons. Such information will help decision-makers to make the decisions for sustainable management in urban environments.

Various researchers proposed different adaptation and mitigation strategies, such as greening the roof [111], increasing urban green space [112], and use of light color or reflective materials to build concrete surfaces, roads, and roofs [113]. Our results agreed with these studies as a positive and negative correlation with LST of built-up areas and green covers, respectively. Furthermore, we need to consider landscape configurations, such as the configurations of green covers and impervious surfaces [7]. Our study also agreed with this study because landscape configurations strongly explained the variability in LST for the different seasons. The urban planner also need to think about the topography and pollution to minimize the urban heat effect as topographic and pollutant parameters significantly affect LST in different seasons.