Spatiotemporal variability of urban heat island: Influence of urbanization on seasonal pattern of land surface temperature in the Metropolitan Region of Belém, Brazil

Spatiotemporal variability of urban heat island: Influence of urbanization on seasonal pattern of land surface temperature in the Metropolitan Region of Belém, Abstract Cities experience the extensive urban heat island effect (UHI), which continue to pose challenges for humanity's increasingly urban population, where tropical cities have experienced a continued and rapid urbanization process in the past few decades. We present the evolution of surface UHI and its controlling factors in the Metropolitan Region of Belém, over the last 16 years (2003 – 2018), which has experienced unique consolidated economic growth and urban transformation under wet equatorial climate. We incorporate MODIS and Landsat satellite data and evaluate statistical techniques for estimates the variation in the land surface temperature (LST) during two seasons: wet season and dry season. Our result revealed that the regions of fast urbanization resulted in a decrease of normalized difference vegetation index and increase of LST. In addition, annual maps showed the spatial pattern of surface UHI intensities were produced based on daytime and nighttime temperature, and the analysis result indicated that the spatial distribution of high heat capacity was closely related with the densely built-up areas. These findings are helpful for understanding the urbanization process as well as urban ecology, which both have significant implications for urban planning and minimize the potential environmental impacts of urbanization in Metropolitan Region of Belém. urbana. Temperatura da superfície terrestre. remoto. Amazônia.


Introduction
Urbanization leads to a dramatic change in the underlying surface structure, properties, and spatial distribution of a city, such as a reduction in green areas and a corresponding increase in impervious areas (Cui et al., 2016;Lee et al., 2019). These changes increase the temperature difference between urban and surrounding non-urbanized areas, or more specifically, indicates that an urban area is significantly warmer than its surrounding rural areas due to artificial land cover and anthropogenic heat (Chakraborty and Lee, 2019), which is known as the Urban Heat Island (UHI). Furthermore, UHI is an important issue for urban planning and environment improvement as it has several implications for energy demand (2), climate adaptation policies (3), public health (4), and heat-related mortality (5) (Paravantis et al., 2017;Agathangelidis et al., 2019). In its broad sense, UHI that is proportional to the degree of urbanization and is closely related to urban climatology, thermal environment, and the quality of human life (Chen et al., 2017;Vahmani et al., 2019). In order to have a comprehensive understanding of the process of urbanization and to evaluate its environmental influence, it is necessary and indispensable to monitor and analyze the dynamics of urbanization in countries.
Traditionally, to address these problems, meteorological departments in many countries regularly measure air temperature at point locations (Shreevastava et al., 2019). However, the complex spatial arrangement of surfaces in urban areas, makes it difficult, or even impossible, to estimate the local variations in surface temperature based on these data alone (Hu et al., 2019). In this way, with the development of remote sensing technology, satellite imagery have been used to estimate Land Surface Temperature (LST), the main driver of air temperature, and to map the spatial distribution of LST, although many studies have been showed the air temperature can be higher or lower than surface temperature depending on various factors such as the presence and direction of wind, insolation and surface characteristics (Lehoczky et al., 2017;Zhou et al., 2018). Nevertheless, LST can provide an estimate of the spatial pattern of temperature over large areas (Mathew et al., 2016).
Thus, the UHI is observed from air temperature measurements within the urban area and the surrounding rural areas based on LST from remote sensors because of high spatial and temporal resolution, free availability, and easy access (Bala et al., 2019). Surface Urban Heat Island (SUHI), namely remotely sensed urban heat island, is usually visualized as a dome of hot air over urban areas with the help of far infrared data that allow to retrieve LST, which they "see" the spatial patterns of upwelling thermal radiance received by the remote sensors (Estoque & Murayama, 2017). This way, LST is the important parameter to analyze urban climatology and is the conventional method to derive the SUHI (Takebayashi and Senoo, 2018;Levermore et al., 2018). It has a direct effect on the air temperature and mean radiant temperature, which relate directly to the thermal environment, and integrating satellite imagery offers numerous opportunities to diagnose potential contributions of physical landscape features that create SUHI (Zhou et al., 2018).
For instance, the Moderate Resolution Imaging Spectroradiometer (MODIS) on board Terra and Aqua, which provides daily LST data with high temporal (four times per day) and relatively high spatial resolution (1 km), has been widely applied for monitoring UHI Niu et al., 2020). Despite having a coarse resolution, the daily MODIS LST products are ideal for time series analysis, providing daily LST data at global scale since 2000. On the other hand, with a high spatial resolution of 30-120 m, Landsat TM/ETM+/OLI data has been used to estimate the LST for SUHI study in recent years (Neiva et al., 2017;Peres et al., 2018;Simwanda et al., 2019). In addition, as vegetation transpiration mitigates the effect of SUHI (Chen et al., 2020), numerous studies have focused on understanding the relationship between LST and the normalized difference vegetation index (NDVI) (Inostroza et al., 2016;Guha and Govil, 2020). These studies found that the NDVI-LST relationship was more suitable for the analysis of SUHI for different seasons in different climate regions (Chakraborty and Lee, 2019). It is therefore increasingly important that long-term data is necessary to capture urban transformation cycles and its relationship with urban development (Hong et al., 2019). In the present study, Landsat 5 TM and Landsat 8 OLI/TIRS thermal images were used to retrieve LST and analyze the spatial pattern of LST and its variation with different LULC types and MODIS image were used to extract SUHI intensity in Metropolitan Region of Belém. Particularly, our analysis focuses on temporal and spatial variability of SUHI intensity along urban development trajectory using MODIS archive and climate conditions (e.g., rainfall) together. The objectives of this study are: (i) to analyze the spatial distributions of NDVI and LST using Landsat TM/OLI/TIRS images acquired in 1988 and 2018; and (ii) to analyze SUHI intensity (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) in Metropolitan Region of Belém.

Study area
Metropolitan Region of Belém is divided into 26 administrative districts, which reflect distinct geographical patterns of climate and weather (IBGE, 2018). The study was conducted in Belém, Guamá, Sacramenta, Entroncamento, Bengui, Icoaraci, Ananindeua, and Marituba districts, which is located in Pará state (01° 30' 00'' to 01° 12' 00'' S latitude and 48° 32' 00'' to 48° 16' 00'' W longitude) (Figure 1), Brazil. It has a total area of 350.34 km², from which the terrestrial part was about 293.35 km². Around the Metropolitan Region of Belém, high-density residential area has been a typical land cover, but the density of residential buildings has been changed dramatically during the study period: compact mid-rise until 1988, and compact low-rise between 1998 and 2018. Therefore, the observed data at the Metropolitan Region of Belém is useful for evaluating and tracking the changes in microclimatic environment with progress of urbanization and economic growth.
The microclimate of area also influences the urban climate (Gonçalves et al., 2019). Tropical cities, and therefore Metropolitan Region of Belém being located in the North Brazil, belong to the tropical wet equatorial zone; af type of Köppen's climate zones (Dubreuil et al., 2018). The elevation ranges between 10 and 40 m.a.s.l., with daily average temperatures around 31.2 and 33.3 °C, and annual precipitation between 2750 and 2834 mm. The bimodal rainy pattern comprises a wet season (December -May) and a dry season (June -November) (Ferreira et al., 2015).

Retrieving land surface temperature
The LST was computed from Landsat 5 TM thermal Band 6, and Landsat 8 OLI/TIRS Band 10 images. To compute LST from these band following four necessary steps are required.

Conversion of Digital Numbers (DN) to spectral radiance
The TM images DN were converted into thermal spectral (Equation 1,2) according to (Markham and Barker, 1986).
where is sensor spectral radiance ( ( 2 * * ) ⁄ ), is quantized calibrated pixel values, is spectral radiance scaled to , is spectral radiance scaled to , is the minimum quantized calibrated pixel values (corresponding to in DN which is 1), is the maximum quantized calibrated pixels values (corresponding to in DN which is 255).
However, DN of Landsat 8 OLI/TIRS of band 10 images were converted to spectral radiance which is given by (Vermote et al., 2016).
where is sensor spectral radiance, is radiance multiplicative scaling factor for the band, is radiance additive scaling factor for the band, is level-pixel value in DN.

Estimation of land surface emissivity
Emissivity refers to the ratio of radiance emitted by real-surface materials at their temperature to radiance emitted by black body at the same temperature. In satellite images, pixels represent the earth surfaces are frequently mixed pixels, that is, they are a mixture of surface types like water, vegetation and soil. To estimate emissivity from satellite thermal band data, NDVI threshold method was used (Wang et al., 2014) by quantifying all land surface types.
Therefore, NDVI and Proportion of vegetation ( ) are required to compute emissivity. The NDVI is an index which infers general conditions of vegetation and used to defined LST. It is also important to compute the proportion of vegetation and emissivity. This index was computed using equation 3.
where is near infrared band reflectance (Band 4 in Landsat 5; and Band 5 in Landsat 8); is red band reflectance (Band 3 in Landsat 5; and Band 4 in Landsat 8).
where is a fully vegetated land-covers NDVI value (0.5); is a non-vegetated landcovers NDVI value (0.2). Then land surface emissivity computed as follows (Sobrino et al., 2004) where is vegetation emissivity; is soil emissivity; is the surface geometrical distribution. Emissivity values of vegetation, soil and water for Landsat 5 and 8 are presented in Table 2.

Computation of brightness temperature
Brightness temperature ( ) is the actual temperature observed by the satellite under emissivity theory (Chander et al., 2009). It can be computed according to Planck's equation (Tan et al., 2009) given below (equation 6). 6/17 where is the effective at sensor brightness temperature in Kelvin; is the spectral radiance at sensor's aperture ( ( 2 * * ) ⁄ ); 1 and 2 are the calibration constant for Landsat; 1 is 607.76 ( 2 * * ) ⁄ ; and 2 is 1260.56 K.

Estimation of LST
The obtained brightness temperatures were referenced to a black body, which differ from the properties of real objects. Correction for spectral emissivity ( ) is a must and the emissivity corrected surface temperature was computed as follows (equation 7): where is the surface radiant temperature in Kelvin (K); is the effective at sensor brightness temperature in Kelvin; is the wavelength of emitted radiance; ℎ is Plank's constant (6.626 * 10 −34 ); is the velocity of light (2.998 * 10 8 ⁄ ); is the Boltzmann constant (1.38 * 10 −23 ⁄ ); and is the surface emissivity.
As the surface radiant temperature is in Kelvin, which differs from the commonly used centigrade, the LST in Celsius was calculated by adding the absolute zero (approximately -273.15 °).

Seasonal variations of SUHI intensity
In this study, to examine the SUHI regarding the spatiotemporal trends of intensity, we defined the urban main built-up area and surrounding non-urbanized area. We then infer built area based on the land cover classification using MapBiomas data (MapBiomas, 2020) and Google Earth Engine (Gorelick et al., 2017). In view of the concentrated distribution of built-up areas in Metropolitan Region of Belém, we extracted the urban main built-up area based on visual interpretation. Water bodies in the urban area were excluded from this analysis because these pixels may have overshadowed the effect of the urbanization on LST (Muro et al., 2018).
We excluded the LST products during the two seasons: May-April and September-November, having very few cloud-free images. Therefore, within each year LST map during the wet season (December -February) and the dry season season (June -August) were considered during 2003 -2018. Specifically, for each year, all the images were divided into daytime (10:00 and 14:00) and night-time (22:00 and 02:00). Then, they were averaged producing two mean LST images in the Metropolitan Region of Belém. In this way, a mean LST pattern for each specific year is available, with nighttime and daytime considered separately, at a 1 -km spatial resolution (Mohammad et al., 2019).
Then, a two-step method was used to estimate the SUHI: (i) the mean daytime and nighttime SUHI intensity, the seasonal and yearly differences in SUHI intensity, the annual maximum and minimum values, their occurring seasons, and the annual ranges of th e SUHI intensity during the daytime, nighttime, and all-day were calculated in the Metropolitan Region of Belém; (ii) significant differences in the SUHI intensity were observed between the daytime and nighttime as well as among the seasons, and these differences were determined using nonparametric tests of two samples, dependent samples, and independent samples. This way, SUHI was defined as the area-averaged differential in LST between urban and rural area (Wu et al., 2019;Lee et al., 2019) as measured by the MODIS-derived LST data from Terra and Aqua satellites as follows: where ∆ represents the difference between the LST in urban and rural area, represents the LST in urban main built-up area, and represents the LST in the rural area (i.e., the background LST). Among them, ∆ 0 indicates the SUHI in the urban area. A positive value of ∆ 0 indicates a SUHI effect, and the opposite indicates a "canopy layer cool island" (SUHI) effect. Therefore, the SUHI is loosely defined as the land surface temperature difference between the urban and rural areas (i.e., ∆ 0 ). Figure 2 shows the derived Landsat TM based composite images, NDVI, and LST. Four reflective multispectral bands (Bands 2-5) were used for compositing. We obtained the spatial pattern of the building coverage in 1988 and in 2018, an obvious urban growth appeared in surroundings non-urban areas in both 1988 to 2018. The value of the NDVI ranges from 0.23 -0.82 in 1988, and changes to 0.13 -0.79 in 2018, and the NDVI changed significantly due to urban expansion during this period. From the spatial distribution in both years, we concluded that the inner-southside area had a smaller NDVI value and a larger value of LST wide-ranging from 28.5 ° -32.7 °, and that through a visual inspection, a surface urban heat island existed within the city.

Spatiotemporal variations of NDVI and LST distribution
The LST spatial pattern shows in 1988 and 2018 considerable differences along the Southeast, North, and West areas of the city. The highest LST values are found in the South -east area of the city, where vegetation cover levels are lower. The lowest values are found towards the northwestern region of the urban area. Lower LST values in the central area of the city may be attributable to higher building densities. Large open spaces have a strong cooling effect. In some cases, cooling islands are associated with higher shares of vegetation, i.e., the presence of parks or other green areas. In other cases, the actual sealed surfaces accounts for lower LST values, probably related to the mutual shadowing of high-rise buildings. It is possible to appreciate the spatial correlation between large open spaces, green areas and lower LST values. A significant correlation was found comparing pixels of LST and NDVI (rPearson = -0.339, p<0.01). The maximum difference in LST based on the coolest and warmest places within the continuous urban area is 8.28 °.

Temporal climate variables and SUHI trends in Metropolitan Region of Belém
Meteorological data from the study period (Figure 3a) obtained from the National Institute of Meteorology (INMET) shows a wide variation in average air temperatures from below to above normal. There was a great variation in annual rainfall during 16 years from the minimum annual rainfall of 2689 mm to the maximum of 3583 mm, and air temperature has shown that its increase is due, at least in part, to differential changes by cloudy conditions, where minimum and maximum air temperature are respectively 31.7 ° and 32.9 °.
In additional, the statistics of diurnal cycle of LST as derived from the average images revealed a strong heat island in the range of 0.85 ° -1.08 ° during day and a weaker heat island in the range of 0.55 ° -1.17 ° at night (Figure 3b). Notably, it is also observed that the rapid increase of SUHI has decreased significantly with a few substantial drops since the late 2010's. Our finding suggests that the urbanization has been moving towards a different phase with completion of urban structure. 9/17

Spatiotemporal pattern of the SUHI in city center
Based on the above variables, the SUHI intensity provides spatially explicit descriptions of the distribution of urban heat throughout the city and by each of the time periods (Figure 4). Land surface temperatures of the output raster SUHI representation depicts average daytime and nighttime maximum surface temperature from 4.23 ° and 1.25 ° with a mean of 0.93 ° and 0.38 °, and standard deviation of 0.84 and 0.44, respectively. In its broad sense, heavily forested areas (including major parks, in addition to certain residential areas) show a tendency towards cooler temperatures. Areas with lower canopy cover, such as those at the northeastern edge of the city (industrial, and airport areas), southeastern freeway corridors, and port yards, consistently appear hotter in the daytime. Additionally, areas around roadways show small and localized 'pockets' of heat within close proximity in nighttime.
We observed a pattern of heat distribution wherein downtown Metropolitan Region of Belém, along with the inner-southside industrial area exhibit the highest levels of heat. Land surface temperatures in these areas can be over 5 ° hotter than areas of the city such as those to the South and North. Therefore, the urban main built-up area shows higher land surface temperatures as compared to the rural reference. Moreover, there is a seasonality influence on the ranking in magnitude of the SUHI intensity between the monitoring stations. In all seasons (wet and dry), however, the built-up locations over South show the highest SUHI intensity values. The differences are larger for the minimum and maximum SUHI values (1.25 ° and 4.35 °) for the mutual differences between the stations. During a large part of the year, SUHI values in the built areas can be substantial ( Figure 5). High values are usually observed on clear (dry season) days. Only in wet season the SUHI intensities are generally lower. However, even in this period SUHI max can reach 4.97 ° or more. While a fast cooling at the rural and water border reference site occurs under these conditions, surface temperatures in the urban environment remain more or less constant. Moreover, in this period, the SUHI intensity reaches its maximum already, remains at a constant maximum during a large part of the daytime, and declines rapidly just before nighttime.
The results clearly reveal that surface temperatures over the south downtown area of Metropolitan Region of Belém were significantly greater than those over the suburbs (North-east area) of Metropolitan Region of Belém in later years. However, Figure 5 reveal that surface temperatures over the Northwest area of Metropolitan Region of Belém are significantly lower than those over the Northeast of Metropolitan Region of Belém on 2005, 2014, and 2017, demonstrating the existence of an "urban cool island" effect, whereby city air is cooler than that of surrounding non-urban areas on daytime. On the other hand, SUHI intensity in Metropolitan Region of Belém tend to be lower in nighttime than in daytime, seems to be similar between the wet and dry as well as on annual (Figure 6). The general trend shows that in all urbanized areas, SUHI intensity distribution tends to be more homogeneous in nighttime, i.e. with a lower land surface temperature amplitude. Unlike, a few areas are recurrent hotspots and coldspots relatively to the spatial mean LST. These climatic inhomogeneities should be related to peculiarities in terms of urban features. In general, due to variations in the SUHI effect and its moderate magnitude during the night, with homogeneous urban morphology and indicators (e.g. inner-southside area) have less coldspots than SUHI intensities with heterogeneous urban morphology (e.g. northwest and northeast areas). Then, clearly that the highest SUHI intensity is distributed in the south downtown and northeast areas, and that lower SUHI intensities are distributed in the northwest and northeast areas between nighttime and daytime in Metropolitan Region of Belém.

Discussion
Accurately patterns of SUHI demands knowledge of how natural and human-created landscapes change weather and climate across spatial and temporal scales (Voelkel and Shandas, 2017;Bauer, 2020). The streetscape characteristics associated with increased SUHI-related daytime and nighttime temperature in Metropolitan Region of Belém were linked to both natural and built causes (Chun and Guldmann, 2018). Across the techniques used in this study, urban features were found to be drastically more important than natural features at higher ambient temperature (Smoliak et al., 2015).
Across study area, common patterns were observed: forested and otherwise vegetated areas are cooler than urbanized areas; lower-density urban areas are often cooler than high-density urban areas. In nighttime, the difference is not changed significantly, in daytime there is the difference is largest, which we speculate may be influenced due to the non-land use variables (e.g., wind speed, albedo, etc.) that are left out of our results (Haashemi et al., 2016). As shown Li et al. (2020), the SUHI phenomenon using satellite remote sensing data have been include variables such as albedo in their techniques which, if included in our research, could improve our analysis power, because many processes observed in nature are spatially variable, especially in the field of geography and meteorology.
In fact, a comparative study of the long-term SUHI phenomenon is not simple because of the complexities regarding acquisition and inconsistency of data. This study, thus, is meaningful because SUHI intensities were compared and analyzed for multiple administrative districts and seasons using unified and available data and methods. MODIS land surface temperature data seem to be able to capture the average increase in temperature as urban areas increase. However, there may be limitations, such as microscopic differences due to the resolution of the data and the use of the temperatures of short-term temperature time series, as demonstrated by Wu et al. (2019).
In additional, De Oliveira et al. (2020) showed that the most negatively evaluated places (e.g., with high temperatures) associate a clear perception of environmental nuisances (e.g., pollution) with a feeling of a degraded or unpleasant condition of the place (e.g. lack of green spaces, natural areas) in Metropolitan Region of Belém. This distinction between places is persistent with season and time of day and is reflected in the evaluation of thermal comfort and consequently lead to increasing effect of SUHI (Priyankara et al., 2019). Thus, loss of urban green areas potentially poses a significant threat to citizens vulnerable to heat-associated illness, particularly considering climate change and global warming predictions, which are set to be exaggerated in a densely built-up area (Venter et al., 2020).
In this way, the stronger negative correlation between mean NDVI and mean LST was also expected, since higher NDVI values signify healthier vegetation, improved photosynthesis, and fuller leaf canopies (Kaplan et al., 2018). In fact, important research remains on how to maximize NDVI values in urban ecosystems, but the literature suggests increasing native species (e.g., flora and fauna) richness, complexity, and diversity (Johnson et al., 2019). However, landscape compositions lack detailed information needed by engineers, architects, urban planners, and designers to create solutions in complex urban settings (Lemonsu et al., 2015;Soltani and Sharifi, 2017).
Based on the correlation analyses, global studies have looked at a limited number of spectral indices, to analyze the influence of land cover on temperature since vegetation cover, surface soil water content and impervious surface cover are known to be determinants of LST (Leconte et al., 2015;Van Hove et al., 2015;Gusso et al., 2017;Alexander, 2020). Although vegetation cover is known to reduce LST, vegetation types may differ in their ability to reduce surface temperature (Sussman et al., 2019) because trees lower surface and air temperature by shading in addition to evapotranspiration. In a study by Duncan et al. (2018), the surface temperature of grass was found to be much lower than asphalt or concrete in sun as well as in shade. Moreover, tree shade had the potential to reduce the temperature of both surfaces by up to 7 °C.
Further studies on urbanization factors and levels influencing the SUHI intensity using higher resolution materials and methods will be needed. Nevertheless, an important yet often unexplored research topic is the spatial patterns of future urban development and its impacts.

Conclusions
Seasonal variation of SUHI were calculated for Metropolitan Region of Belém for the period 2003-2018, using MODIS and Landsat satellite data. LST and land cover data were used in order to assess the LST values in urban and non-urban areas of the city and to examine the SUHI intensity trend for the last 16 years. The 2003-2018 time-series analysis found positive daytime LST trends, with the highest values found at the urban core of the city, and positive nighttime LST trends. SUHI trends exhibit large variations, for daytime an increasing but not statistically significant. At nighttime, the SUHI trends are considerably limited with both positive and negative trends found.
Our results highlight the important contribution of remote sensing towards improving the maps of urban heat islands in tropical areas. Using satellite imagery to expand insights for SUHI is particularly useful in data-poor areas and should be helpful for policymakers to formulate countermeasures to mitigate the effects of SUHI and create more sustainable and environment-friendly cities. In addition, further comparison studies among different cities could be conducted in order to obtain more reliable results, and different climatic conditions should also be considered in different cities in future research.