Spatiotemporal Analysis of Urban Heat Islands and Vegetation Cover Using Emerging Hotspot Analysis in a Humid Subtropical Climate

: Research on the temporal and spatial changes of the urban heat island effect can help us better understand how urbanization, climate change, and the environment are interconnected. This study uses a spatiotemporal analysis method that couples the Emerging Hot Spot Analysis (EHSA) technique with the Mann–Kendall technique. The method is applied to determine the intensity of the heat island effect in humid subtropical climates over time and space. The data used in this research include thermal bands, red band (RED) and near-infrared band (NIR), and Landsat 7 and 8 satellites, which were selected from 2000 to 2022 for the city of Sari, an Iranian city on the Caspian Sea. Pre-processed spectral bands from the ‘Google Earth Engine’ database were used to estimate the land surface temperature. The land surface temperature difference between the urban environment and the outer buffer (1500 m) was modeled and simulated. The results of this paper show the accuracy and novelty of using Emerging Hotspot Analysis to evaluate the effect of vegetation cover on the urban heat island intensity. Based on the Normalized Difference Vegetation Index (NDVI), the city’s land surface temperature increased by approximately 0.30 ◦ C between 2011 and 2022 compared to 2001 to 2010. However, the intensity of the urban heat island decreased during the study period, with r = − 0.42, so an average − 0.031 ◦ C/decade decrease has been experienced. The methodology can be transferred to other cities to evaluate the role of urban green spaces in reducing heat stress and to estimate the heat budget based on historical observations.


Introduction
For several reasons, understanding temporal and spatial variations in urban heat islands (UHI) is critical [1].First, the UHI effect is a phenomenon in which metropolitan regions have greater temperatures than surrounding rural areas due to human activities such as transportation, industrial processes, and building materials [2].This effect can negatively influence human health (heat-related illnesses and mortality) and the environment (increased energy usage and greenhouse gas emissions) [3].Understanding the temporal and spatial variations of the heat island effect can thus assist policymakers and urban planners in developing effective mitigation solutions to mitigate its detrimental effects [4].
Second, the UHI effect can have economic consequences, such as increased energy costs for cooling buildings and damage to infrastructure due to extreme heat [5].The UHI effect is a complex phenomenon influenced by local and global influences.Green areas and water-consuming urban heat island effect (UHIE) solutions can help mitigate the consequences of the UHI effect.Still, more research is needed to fully understand the influence of global warming on the UHI effect.However, urban activities such as transportation, buildings, industry, and deforestation contribute to climate change, and measuring climate change awareness in the public is critical for designing mitigation measures [6].In general, the findings of various studies on the intensity of urban heat islands in various dimensions show that urban heat islands are a phenomenon in which urban areas tend to be warmer than their rural counterparts due to heat absorption and retention by buildings and other infrastructure [7].
Third, researching the temporal and spatial variations of the UHI effect can benefit from new spatiotemporal analysis methods and available data of geospatial information based on satellite maps and remote sensing.Urban planners and policymakers can benefit from the advancements of Geographical Information Systems (GIS) and earth observation satellites to identify locations most sensitive to these economic repercussions and design targeted mitigation plans by examining the temporal and spatial fluctuations of the UHI effect.Research on UHI effects using satellite map datasets with historical observations of land surface temperatures, thermal bands, red bands, and near-infrared bands is proliferating [8][9][10].Emerging Hot Spot Analysis (EHSA) is a technique that falls under Exploratory Spatial Data Analysis (ESDA) hot spot analysis.The Emerging Hot Spot Analysis tool calculates the Getis-Ord Gi* statistic using the traditional time-series Mann-Kendall test for monotonic trends.
Therefore, this study uses a spatiotemporal analysis method to determine the intensity of the heat island effect in humid subtropical climates over time and space.To evaluate whether urbanization and population expansion or other reasons, such as climate change, are to blame for temporal changes in UHII (Urban Heat Island intensity), it is crucial first to comprehend the changes in UHI's temporal and spatial pattern using a case study.The case study was chosen as a mid-size study representative of a city in a humid subtropical climate (Köppen: Cfa, Trewartha: Cf).Sari, located in Northern Iran, was chosen as a representative city on the Caspian Sea with high population growth rates.Between 1960 and 2020, the city population increased 17 times, reaching 338 thousand people in 2025 [11].While the impact of climate change on the alteration of air temperature patterns in the city of Sari has been assessed previously, a review of the existing literature reveals a gap in research regarding the influence of global warming on changes in the Urban Heat Island Intensity (UHII) for this city.Consequently, the discoveries and framework presented in this study can serve as a valuable foundation for future articles aiming to investigate and assess the impact of global warming on coastal cities in northern Iran.The research questions of this research are formulated as follows:

•
How to couple Emerging Hot Spot Analysis and the time-series Mann-Kendall test for characterizing the vegetation change and urban heat island effects?

•
What are the effects of climate change on the urban heat island effect intensity in Sari?

•
What is the influence of vegetation on the land-surface temperature in Sari?

Literature Review
A UHI is where urban regions suffer higher temperatures than surrounding rural areas due to human activities and changes in urban composition and layout [12].Due to human activities such as population pressure, infrastructure development, and economic growth, urban regions experience greater temperatures than neighboring rural areas [13].UHI has the strongest effect on large and densely populated cities and impacts many human societies worldwide [7,14].The trend of urban integration has accelerated, amplifying the impacts of UHI.UHI impacts are caused by the interaction of land use pattern change, subsurface diversity, and human activities [15].The energy foundation of the urban heat island can be linked to gravity, which is extremely powerful since matter is concentrated in a relatively small region [7,16].The UHI effect is a severe environmental and social issue that modern cities worldwide face, posing threats to urban populations' physical and mental health, air quality, and the healthy growth of surface plants.This issue will grow more pressing as urbanization continues [17].UHI is caused by urbanization and leads to changes in land use, vegetation degradation, and climate change [13].
Global warming can impact UHI since an increase in urban air temperature due to UHI or global warming can be advantageous in the winter but detrimental in the summer [16].The findings of a study on the influence of global warming on the UHI demonstrate that in some circumstances, modern climate change can increase and decrease the severity of the UHI effect [18].Urban stations frequently demonstrate local UHI effects and global warming, making it difficult to distinguish between the two signals in the reported surface air temperature (SAT) warming.The mean urban impact indicator (Uii) of the 45 urban stations increased from 0.06 to 0.35 from 1992-2013; this increase was revealed by linear regressions between the trends of annual averaged daily mean SAT, maximum SAT, and minimum SAT between urban and near stations [19].Additionally, positive correlations between the changes in Uii and SAT trends were substantial at the 0.001 significance level.The combination of heat waves (HWs) and UHI has increased the intensity, frequency, and duration of extreme heat events in urban areas, seriously endangering the health of urban populations.According to Luo et al. (2023) [20], low wind speeds contribute to increased HW frequency of HWs and augmentation of UHII.Both the UHII and HWs decreased with rising wind speed.Anthropogenic heat (AH) emissions from human activities can alter urban circulation, which in turn affects the level of air pollution in and around cities. AH emissions can considerably alter urban heat islands and urban breeze circulations in large cities.These variations in weather can considerably impact the geographical and vertical distributions of air contaminants [21].Urbanization and economic growth are responsible for the deterioration of the environment and lead to increased carbon emissions [22].Urbanization causes thermal pollution caused by specific urban activities and patterns on the ground levels in cities, which changes the city's local climate.It at least negatively affects the comfort of the city [23].Based on investigations in six major cities surrounding the Hangzhou Gulf in eastern China, the climatic effects of the heat island's growth and intensity suggest that the phenomenon is increasing irregularly in the research area.
Because solar and human activities impact UHI days, the change during the day appears bigger than the change at night.There was some urban cool island phenomenon in winter when urban surface temperatures were −1-0 • C lower than the ambient LST.Except for winter, other seasons' UHI during daytime changed from 0.5 • C to 2.5 • C and during nighttime strongly changed from 0.5 • C to 5 • C in 2009 [24].A preliminary study of land temperature in the Bangkok metropolitan area revealed a trend toward higher ground temperatures in several seasonal sectors, especially agricultural activities.Burning straw and weed cover during the preparation of rice cultivation land caused the greatest LST in the rice fields [25].According to a study by Nyangena et al. (2019) [26], sub-Saharan Africa's rapid urbanization and economic growth resulted in higher demands for housing, transportation, and energy.According to the study, policies, and strategies considering environmental issues and carbon emissions are required to support sustainable urbanization and economic growth.Finally, the authors concluded that urbanization has a big impact on climate change and that urbanization and how nations may develop without affecting the environment should be the focus of future research.To determine how urban and physical characteristics affect air temperatures and produce different UHI configurations, an experimental study was conducted at the micro-scale at several spots in the city.The intensity of UHI peaks during the summer and at night.Finally, an increase of up to 4 • C in urban air temperatures will be expected, significantly increasing urban livability, according to Martinelli et al., 2020 [27], who studied the frequency distribution of current and future weather situations.UHII in Coventry showed that UHI intensity does not follow a predictable pattern across the seasons.However, minor temperature rises did occur in some urban locations [28].The impact of UHI on the heating intensity of buildings was analyzed based on the observed hourly temperature data during the heating periods from 2009 to 2017 in Tianjin, China.The results showed that the average annual heating intensity decreased by 8.08% and 10.30%, with the urban center's temperature increasing by 1.46 • C and 1.90 • C compared with the suburban and outer suburban areas [29].The rise in temperatures within the Mediterranean region can worsen the UHI.
Consequently, cities like Rome and Thessaloniki experience minimum temperatures (Tmin) compared to rural areas, where nighttime UHIs range from +1.5-3 • C.However, researchers did not find any changes in UHI due to the impact of climate change [30].In the most disastrous climate change scenarios, the number of heat wave days might treble if there is no UHI effect, according to a study carried out for Sydney, Australia, because, on more than 90% of the days, the heat wave does not pass the threshold [31].
Therefore, this paper aims to help us better comprehend the complex relationships between urbanization, climate change, and the environment.This knowledge can be used to create sustainable urban planning and design methods that can help alleviate urbanization's negative environmental impacts and promote more sustainable and resilient communities.Overall, research into the temporal and spatial changes of the heat island effect is vital for understanding the effects of urbanization on the environment, human health, and the economy and creating effective mitigation strategies.

Materials and Methods
Sari is situated at a geographical location of 53 • 3 E longitude and 36 • 34 N latitude.It has an average elevation of approximately 40 m above sea level.Sari's general topography falls within the 0-100 m altitude range.The northernmost part of the city stands at a height of 12 m, while the southernmost part reaches 80 m.The alignment lines in Sari run from east to west.
The city experiences a warm and humid subtropical climate (Köppen: Cfa, Trewartha: Cf), which shares similarities with a Mediterranean climate (Csa).Winters in Sari are cool and rainy, while summers are hot and humid.With an average annual temperature of 18 • C, it is slightly warmer than the Iranian average.While July and August bring scorching heat with highs above 30 • C, January and February offer a respite with lows around 4 • C. Rainfall averages around 665 mm annually, peaking in spring and fall, although recent trends indicate a slight decline.High humidity, particularly in summer (up to 80%), adds to the perceived warmth.Moderate winds shift direction with the seasons, blowing from the northwest in winter and the southwest in summer.Sari has over 2183 h of sunshine annually, which is more than in other parts of Iran [32,33].
Sari has a rich history that dates to the Sassanid era and has evolved substantially over time [34].Due to fast population expansion and urbanization, Sari has expanded its city borders and created new residential and commercial areas in recent decades (Figure 1).Sari's strategic location as a transportation hub in northern Iran has played a significant role in its growth.This city has a railway station and an airport, and it is located on the main highway that connects Tehran to the Caspian Sea.Another factor impacting the development is that Sari is a hub for regional businesses, health services, and education.The city's numerous educational institutions, medical facilities, and commercial centers attract scholars, patients, and customers from the neighboring districts [35,36].
In this study, Landsat satellites 7 and 8 from 2000 to 2022 collected thermal bands, red bands (RED), and near-infrared bands (NIR).For this purpose, we considered days with 1% or less cloud cover (Table 1).
The LST was estimated using pre-processed spectral bands in a 'Google Earth Engine' environment.The BT parameter (brightness temperature) was determined by analyzing the thermal band.Land Surface Emissivity (LSE) was then estimated using the NDVI index.For this purpose, NDVI was first estimated using Equation (1) [37,38]: In this study, Landsat satellites 7 and 8 from 2000 to 2022 collected thermal bands, red bands (RED), and near-infrared bands (NIR).For this purpose, we considered days with 1% or less cloud cover (Table 1).In Table 2, the details of urban vegetation classes and NDVI value are provided [39]: Equation ( 2) was used to calculate the proportion of vegetation (PV) [40].
In Equation ( 4), W is equivalent to the wavelength of emitted radiance.The LST difference between the urban area and the surroundings was measured for the UHII of Sari City in the following section following the LST assessment in the study region.This study chose a 1500 m border in Sari's outer environment (buffer) as a rural area.Based on this, the following equation [41] was also extracted to produce the UHII: In Equation (5), Turban stands for the LST in urban areas, and Trural stands for the LST in the city's outside regions.Then, using ArcGisPro software (Version 10.8.2), the needed data in the form of structure was used to identify the trend of changes and the zoning of the anomalies of UHII over Sari.It should be mentioned that the Mann-Kendall technique was utilized for the temporal and spatial analysis of the trend of changes in the UHII of the city of Sari, and the hotspot analysis function was employed for the zoning of the intensity anomaly of the UHI.This method calculates the difference between each observation and all the observations, and the parameter S is obtained.For a random sample containing n observations, the estimator S can be calculated with the following equation [42,43]: where n is the total number of series observations, and x j and x k indicate the jth and kth data points, respectively.The "sgn" function may alternatively be computed using the Equation ( 7) [44][45][46]: Next, one of the following equations (Equation ( 8)) is used to calculate the variance of S [44][45][46]: where n and m indicate how many sequences include at least one duplicate data value.
The frequency of similar data in a sequence is also represented by t.Finally, one of the following relations (Equation ( 9)) is used to obtain the Z [44][45][46]: In the case of a two-tailed trend test, the following criteria must be satisfied to accept the null hypothesis: |Z| < Z a/2 (10) where a is the significance level that is considered for the test and Z a is the standard normal distribution statistic at the significance level of a, which becomes a 2 because of the two domains of the test.
It is necessary to explain that in this study, the main pixels throughout the study area are combined in the form of 3 × 3. Therefore, to perform this process, the fishnet function was used to obtain pixels with dimensions of 88 m.Then, the significance of the change process was evaluated using the Mann-Kendall method for each of the pixels.
The UHII anomaly was then estimated using the Getis-Ord Gi statistic and is determined as follows [41]: In Equation (12), Xj is the attribute value for complication j, Wij is the spatial weight between complications i and j, and n is equal to the total number of complications.
Gi is some Z-score.Therefore, no further estimate is required.Subsequently, we applied the Emerging Hotspot Analysis (EHA) using ArcGIS Pro software [47] to understand spatiotemporal variations of emergency calls during the study period.The Emerging Hotspot Analysis is among the new methods added to the GIS-based analyses, and its usage in spatiotemporal analysis is growing [48][49][50][51][52][53][54].This information is used to analyze various zones according to how frequently EHSA occurred in various metropolitan locations throughout the 23-year timeframe.Using ArcGisPro3.1 software, the classes of this assessment, according to Table 3, were employed.Table 3. Spatial and temporal patterns for Emerging Hot Spot Analysis [47].

Pattern Name Definition
No Pattern Detected Does not fall into any of the hot or cold spot patterns defined below.
New Hot Spot A location that is a statistically significant hot spot for the final time step has never been a statistically significant hot spot.

Consecutive Hot Spot
A location with a single uninterrupted run of statistically significant hot spot bins in the final time-step intervals.The location has never been a statistically significant hot spot before the final hot spot run, and less than ninety percent of all bins are statistically significant hot spots.

Intensifying Hot Spot
A location that has been a statistically significant hot spot for ninety percent of the time-step intervals, including the final time step.In addition, the intensity of clustering of high counts in each time step is increasing overall, which is statistically significant.

Pattern Name Definition
Persistent Hot Spot A location that has been a statistically significant hot spot for ninety percent of the time-step intervals with no discernible trend indicating an increase or decrease in the intensity of clustering over time.

Diminishing Hot Spot
A location that has been a statistically significant hot spot for ninety percent of the time-step intervals, including the final time step.In addition, the intensity of clustering in each time step is decreasing overall, and that decrease is statistically significant.

Sporadic Hot Spot
A location that is an on-again, off-again hot spot.Less than ninety percent of the time-step intervals have been statistically significant hot spots, and none have been statistically significant cold spots.

Oscillating Hot Spot
A statistically significant hot spot for the final time-step interval with a history of also being a significant cold spot during a prior time step.Less than ninety percent of the time-step intervals have been statistically significant hot spots.

Historical Hot Spot
The most recent period is not hot, but at least ninety percent of the time-step intervals have been statistically significant hot spots.
New Cold Spot A location that is a statistically significant cold spot for the final time step and has never been a statistically significant cold spot before.

Consecutive Cold Spot
A location with a single uninterrupted run of statistically significant cold spot bins in the final time-step intervals.The location has never been a statistically significant cold spot before the final cold spot run, and less than ninety percent of all bins are statistically significant cold spots.

Intensifying Cold Spot
A location that has been a statistically significant cold spot for ninety percent of the time-step intervals, including the final time step.In addition, the intensity of clustering of low counts in each time step is increasing overall, which is statistically significant.

Persistent Cold Spot
A location that has been a statistically significant cold spot for ninety percent of the time-step intervals with no discernible trend, indicating an increase or decrease in the intensity of clustering of counts over time.

Diminishing Cold Spot
A location that has been a statistically significant cold spot for ninety percent of the time-step intervals, including the final time step.In addition, the intensity of clustering of low counts in each time step is decreasing overall, and that decrease is statistically significant.

Sporadic Cold Spot
A location that is an on-again, off-again cold spot.Less than ninety percent of the time-step intervals have been statistically significant cold spots, and none have been statistically significant hot spots.

Oscillating Cold Spot
A statistically significant cold spot for the final time-step interval with a history of being a statistically significant hot spot during a prior time step.Less than ninety percent of the time-step intervals have been statistically significant cold spots.

Historical Cold Spot
The most recent period is not cold, but at least ninety percent of the time-step intervals have been statistically significant cold spots.

Analyses of the LST Changes and UHII
Point-to-point LST differences were obtained using the average rural area in this investigation to explore the temporal and spatial changes of the UHI.The UHI in the city was used for demonstrating regions warmer than 0 • C. Figure 2a displays the change in the average LST rural and the area of Sari city during the 2001-2022 period.As can be observed, the LST in the city region was higher than the LST in rural areas over the entire period.However, its alterations are accompanied by a growing tendency and fluctuation in different years.However, Figure 2b depicts the difference between the average yearly LST in rural and urban areas, which shows a considerably significant decreasing trend (r = −0.42).In other words, the UHI's severity has been lowered due to urban warming.
It appears that the escalation of global warming is not confined solely to urban regions but extends to suburban areas as well.Consequently, this phenomenon results in a reduction of temperature disparity between the city and its suburbs, diminishing the intensity of the UHI.Despite this challenge, certain mechanisms contribute to UHI intensity reduction, with one notable factor being alterations in urban green spaces, as expounded upon in the subsequent sections.
A parallel issue is observable when evaluating temperature differences between elevated and low-lying regions.The impact of global warming manifests in a diminished temperature contrast between high-altitude mountainous areas and low-altitude areas [55,56].It appears that the escalation of global warming is not confined solely to urban regions but extends to suburban areas as well.Consequently, this phenomenon results in a reduction of temperature disparity between the city and its suburbs, diminishing the intensity of the UHI.Despite this challenge, certain mechanisms contribute to UHI intensity reduction, with one notable factor being alterations in urban green spaces, as expounded upon in the subsequent sections.
A parallel issue is observable when evaluating temperature differences between elevated and low-lying regions.The impact of global warming manifests in a diminished temperature contrast between high-altitude mountainous areas and low-altitude areas [55,56].°C and 29.81-29.9°C.This range is equivalent to 80.15% of the city's area in the first period and 94.26% of Sari's area in the second period .However, the LST class range increased between 2010 and 2015, and the greatest urban area was between LST classes 29.9-29.81°C and 30.30-30.21°C.So, from 2010 to 2015, the area of coverage of these two tiers covered 97.19% of the city, and from 2016 to 2020, it covered 96.44% of the city.According to Figure 3, the range of LST changes in the city is generally increasing .Figure 5 shows the area changes of the UHII in different classes and the skewness of the UHII trends to be positive in 2001-2005 compared to 2016-2020.In other words, the greater part of city areas with decreased UHII occurred between 2016 and 2020.An important sign of the UHII is that, in addition to existing in areas with low construction density, the UHI has not spread.Rather, the Tajen River's role in weakening the power of the UHI is pretty evident, and it has remained so across all four periods.In general, the results of this phase of the study indicate a decrease in UHII with the rise in NDVI in recent years.These trends are particularly noticeable in the second through fourth decades.The overall average NDVI for these decades is 0.211, 0.228, and 0.238, respectively.Correspondingly, the average UHII for the second and third decades is 0.197 °C, while for the fourth decade, it is 0.172 °C.Conversely, the findings reveal an increasing trend in both high vegetation and low vegetation areas from the second to the fourth period.This growth in green spaces within the urban environment directly contributes to the reduction in UHII.In other words, the greater part of city areas with decreased UHII occurred between 2016 and 2020.An important sign of the UHII is that, in addition to existing in areas with low construction density, the UHI has not spread.Rather, the Tajen River's role in weakening the power of the UHI is pretty evident, and it has remained so across all four periods.In general, the results of this phase of the study indicate a decrease in UHII with the rise in NDVI in recent years.These trends are particularly noticeable in the second through fourth decades.The overall average NDVI for these decades is 0.211, 0.228, and 0.238, respectively.Correspondingly, the average UHII for the second and third decades is 0.197 • C, while for the fourth decade, it is 0.172 • C. Conversely, the findings reveal an increasing trend in both high vegetation and low vegetation areas from the second to the fourth period.This growth in green spaces within the urban environment directly contributes to the reduction in UHII.

Investigating the Trend of Changes in the UHII of Sari City
Figure 6 shows a distribution map of the significant change in the intensity of Sari's UHI from 2000 to 2022.As previously indicated, the Mann-Kendall approach was utilized to create this estimate as independent pixels in a 23-year time series.As seen on the map, the main parts of the city center indicate a downward trend, while the city's outskirts show an increase in UHII (upward trends).This issue indicates that the city's outskirts are fast developing.Table 4 demonstrates that, on the other hand, about 11.99% of the city area with a 99% confidence level and 5.77% with a 95% confidence level reveal an upward trend of UHII in Sari City.While 27.14% of the city area has a confidence level of 99%, and 8.81% has a confidence level of 95%, the severity of the UHII is decreasing (downward trend).It has a larger overall extent than the expanding state in the decreasing trend areas of the UHI.

Investigating the Trend of Changes in the UHII of Sari City
Figure 6 shows a distribution map of the significant change in the intensity of Sari's UHI from 2000 to 2022.As previously indicated, the Mann-Kendall approach was utilized to create this estimate as independent pixels in a 23-year time series.As seen on the map, the main parts of the city center indicate a downward trend, while the city's outskirts show an increase in UHII (upward trends).This issue indicates that the city's outskirts are fast developing.Table 4 demonstrates that, on the other hand, about 11.99% of the city area with a 99% confidence level and 5.77% with a 95% confidence level reveal an upward trend of UHII in Sari City.While 27.14% of the city area has a confidence level of 99%, and 8.81% has a confidence level of 95%, the severity of the UHII is decreasing (downward trend).It has a larger overall extent than the expanding state in the decreasing trend areas of the UHI.

Investigating the Trend of Changes in the UHII of Sari City
Figure 6 shows a distribution map of the significant change in the intensity of Sari's UHI from 2000 to 2022.As previously indicated, the Mann-Kendall approach was utilized to create this estimate as independent pixels in a 23-year time series.As seen on the map, the main parts of the city center indicate a downward trend, while the city's outskirts show an increase in UHII (upward trends).This issue indicates that the city's outskirts are fast developing.Table 4 demonstrates that, on the other hand, about 11.99% of the city area with a 99% confidence level and 5.77% with a 95% confidence level reveal an upward trend of UHII in Sari City.While 27.14% of the city area has a confidence level of 99%, and 8.81% has a confidence level of 95%, the severity of the UHII is decreasing (downward trend).It has a larger overall extent than the expanding state in the decreasing trend areas of the UHI.

Investigating the Trend of Changes in the NDVI of Sari City
Figure 7 shows a distribution map of the significant change in the NDVI over Sari from 2000 to 2022.As seen on the map, the main parts of the city center indicate an upward trend, while moving away from the city center, the NDVI change trend is decreasing (downward trend).This can be because, with the city's expansion (urban sprawl), the surrounding lands have been allocated to urban use, and on the other hand, the green space in the city's central areas has been expanding.Table 5 demonstrates that, on the other hand, about 12.06% of the city area with a 99% confidence level and 3.60% with a 95% confidence level reveal an upward trend of NDVI in Sari City.While 5.85% of the city area has a confidence level of 99%, and 11.84% has a confidence level of 95%, the severity of the NDVI is decreasing (downward trend).In general, the area with a significant level of more than 90% for the increasing trend of NDVI with an area of 48.6% is more than the areas with a decreasing trend of NDVI with an area of 18.2% for the city of Sari.
the other hand, about 12.06% of the city area with a 99% confidence level and 3.60% with a 95% confidence level reveal an upward trend of NDVI in Sari City.While 5.85% of the city area has a confidence level of 99%, and 11.84% has a confidence level of 95%, the severity of the NDVI is decreasing (downward trend).In general, the area with a significant level of more than 90% for the increasing trend of NDVI with an area of 48.6% is more than the areas with a decreasing trend of NDVI with an area of 18.2% for the city of Sari.Based on the long-term average from 2000 to 2020, 0.22% of the city's area is included in the water body class.The No Vegetation class occupies the largest area with 51.83%, and for the two classes of Low Vegetation and High Vegetation, the areas are 43.82% and 4.13%, respectively (Figures 4 and 8).In general, in every decade, on average, the area of the water body class has increased by 0.23%/decade, which is a very small percentage.But it is interesting to note that the area of the No Vegetation class has decreased by −9.97%/decade; on the other hand, the area of the Low Vegetation class shows an average increase of 10.85%/decade.At the same time, the areas covered by the High Vegetation class, which are usually around the city and are the place of physical expansion and development, show a decrease of −1.11%/decade (Figures 4 and 8).

Temporal and Spatial Changes of Emerging Hot Spots of Urban Heat Island in Sari City
The EHSA tool was applied in this research to assess spatiotemporal patterns of UHII anomalies within Sari.As mentioned earlier, the tool utilizes two statistical measures: first, the Getis-Ord Gi* statistic and second, the Mann-Kendall trend test.The Mann-Kendall statistic is used to determine if there is a statistically significant temporal trend in each bin's 23-year time series of Z scores from the Getis-Ord Gi* statistic.
However, it also reveals numerous Z scores for the complexity over time in hot spot time series analysis.As a result, distinct locations are given different names based on the frequency at which they occur.Figure 9 illustrates the findings of the EHSA analysis function in this investigation over 23 years.Based on this, the hot spot pattern recognized approximately 31.6% of the entire city, the cold spot pattern identified 49.26%, and the lack of pattern identified 19.24%.However, according to Table 6, the "Oscillating Cold Spot" accounts for 21.8% of the city in the cold spot pattern.In addition, "Persistent Cold Spot" extends to 11.6% of the city area in this manner.Based on the long-term average from 2000 to 2020, 0.22% of the city's area is included in the water body class.The No Vegetation class occupies the largest area with 51.83%, and for the two classes of Low Vegetation and High Vegetation, the areas are 43.82% and 4.13%, respectively (Figures 4 and 8).In general, in every decade, on average, the area of the water body class has increased by 0.23%/decade, which is a very small percentage.But it is interesting to note that the area of the No Vegetation class has decreased by −9.97%/decade; on the other hand, the area of the Low Vegetation class shows an average increase of 10.85%/decade.At the same time, the areas covered by the High Vegetation class, which are usually around the city and are the place of physical expansion and development, show a decrease of −1.11%/decade (Figures 4 and 8).

Temporal and Spatial Changes of Emerging Hot Spots of Urban Heat Island in Sari City
The EHSA tool was applied in this research to assess spatiotemporal patterns of UHII anomalies within Sari.As mentioned earlier, the tool utilizes two statistical measures: first, the Getis-Ord Gi* statistic and second, the Mann-Kendall trend test.The Mann-Kendall statistic is used to determine if there is a statistically significant temporal trend in each bin's 23-year time series of Z scores from the Getis-Ord Gi* statistic.
However, it also reveals numerous Z scores for the complexity over time in hot spot time series analysis.As a result, distinct locations are given different names based on the frequency at which they occur.Figure 9 illustrates the findings of the EHSA analysis function in this investigation over 23 years.Based on this, the hot spot pattern recognized approximately 31.6% of the entire city, the cold spot pattern identified 49.26%, and the lack of pattern identified 19.24%.However, according to Table 6, the "Oscillating Cold Spot" accounts for 21.8% of the city in the cold spot pattern.In addition, "Persistent Cold Spot" extends to 11.6% of the city area in this manner.Moreover, the "Intensifying Cold Spot" has extended to 7.29% of the city.Then, it should be mentioned that 5.86% of the city of Sari had the "Diminishing Cold Spot" pattern present.Cold spot occurrence zones with the "Sporadic Cold Spot" pattern can be seen in 3.02% of the city's area.But based on the distribution of hot spots, "Diminishing Hot Spots" occupy 8.73% of the area, and the "Oscillating Hot Spot" pattern is equivalent to 8.54% of the city area.

Discussion
This study aimed to characterize and comprehend the changes in urban heat island temporal and spatial patterns in Sari, Iran, as a representative humid subtropical climate.The study investigated vegetation changes as a factor influencing the behavior of UHII.In the following paragraphs, we discuss the major study findings and limitations and reflect on the requirements for future work.
The intensity of the UHI is one of the climatic phenomena influenced by climate change and can also impact local climates.In other words, the large-scale climate can influence the UHII while controlling the small-scale climate.The average temperature of the metropolitan environment and its outer space are a function of an increasing trend in global temperature, according to the estimated results of this study since both habitats exhibit an increasing temperature trend.In the case of Sari, Iran, the difference between the metropolitan environment's average LST and its outermost boundary is decreasing.In other words, it can be inferred that the temperature rise brought on by global warming decreases the temperature disparity between tiny climate regions.The LST range between 2011 and 2022 has expanded by around 0.30 As a result, it can be thought that the LST of Sari has been increasing in recent years due to global warming.Other variables, such as human indicators (AI), could have influenced it.
On the other hand, the examination of the UHII in Sari suggests that the range of UHI classes is the same in all 5-year periods, except that the class area has been relocated.Thus, the largest area from 2001 to 2005 is 0.30-0.21• C, accounting for 32.2% of the total city, whereas the largest area from 2016 to 2020 is in the range of 0.20-0.11• C, accounting for 32.2%.In other words, even though the total severity of the UHI has not changed in 20 years, its placement inside the urban boundaries has moved.Sari's urban development and extension have been affected by time and place.Climate change can justify LST fluctuations based on LST variations and the UHII.
Uncertainty in remote sensing data, particularly in the context of NDVI data, is a critical consideration that can impact the reliability and interpretation of results.The sources of uncertainty in NDVI data are diverse, ranging from sensor limitations and atmospheric conditions to preprocessing techniques.Variability in sensor calibration, spectral resolution, and radiometric accuracy introduces uncertainties that may propagate throughout the entire data processing chain.Additionally, atmospheric conditions, such as cloud cover, aerosols, and water vapor, can affect the quality of satellite imagery and introduce uncertainties in the derived NDVI values.Furthermore, choices made during data preprocessing, such as atmospheric correction and geometric rectification, can influence the final accuracy of NDVI data.Acknowledging and quantifying these uncertainties are crucial for ensuring the robustness of any analysis or modeling based on remote sensing-derived NDVI data.Researchers often employ sensitivity analyses and validation techniques to understand better and communicate the potential sources and magnitude of uncertainty in their findings [57,58].
Our study evaluated the NDVI variability and monitored it during the last two decades.Landsat satellites 7 and 8 data were from 2000 to 2020.Thus, our research methodology is robust and valid based on a recent dataset of land surface temperatures.Moreover, we applied state-of-the-art hotspot analysis to understand the spatiotemporal variations.Energy Hot Spot analysis is a new method used in GIS-based analysis.The richness and freshness of the dataset allowed us to draw a novel and accurate picture of the UHI in Sari, Iran.The satellite data allowed us to calculate and visualize the average LST difference between the city and the surrounding area.However, the temporal and spatial dispersion of Sari UHI can be influenced by human activity.This conclusion is supported by the fact that the intensity range of the UHI at the surface has not changed.However, its average value shows a decreasing tendency, indicating that it is a function of climate change, which causes LST uniformity.The Mann-Kendall approach reveals a considerable drop in the UHI's severity in Sari's center area.
In point-to-point urban mode, 35.95% of the urban area decreased, although the marginal expansion of the UHII increased in recent years.This research showed that the increase in the UHII in suburban areas is due to decreased NDVI and the decrease in the High Vegetation class.In fact, with the expansion of urban areas to the surrounding areas, urban land use has replaced the dense vegetation cover, and the UHII has increased in this area.On the other hand, the decrease in UHII in the central areas is also affected by the increase in NDVI, especially the increase in the Low Vegetation class.
Spatial and temporal analysis of the urban heat island effect in cities is gaining attention [59].Several studies are reporting the results due to the availability of remotesensing-based datasets [60].The increase in the surface of the urban heat island in many cities, including Sari, is becoming a common finding [61].Most of those studies focus on coupling land surface temperature and NDVI [62].Our results and methodology benefit from the advancement of spatio-temporal modeling utilizing remote sensing techniques and geostatistical hot spot analysis [63].However, what is needed is to find a relationship between the surface heat island effects and air heat island effects [64].Remote sensing can stay on the 2D level if we do not address the air heat temperature and think about the urban thermal layer in 3D [65].
Future research is needed to properly comprehend the connection between the NDVI and UHI impact mitigation.By increasing the green space in the urban area, the intensity of thermal stress and the temperature budget of the city can be reduced.Overall, the impact of introducing green vegetation or nature-based solutions into the UHI effect is complex and varies depending on the city's location, layout, and growth, the quantity and density of its inhabitants, and the urban economy's energy and water use.As shown in this study, increasing green spaces in cities can counterbalance the effects of the UHI.Also, field surveys and local measurement campaigns are needed to move from the mesoscale to the microscale to move to a lower urban thermal layer, addressing the spatial and temporal variability [66].Decisions regarding how to adapt to control urban densification and lessen the effects of climate change in urban settings can be made with the use of this knowledge by policymakers and urban planners.There is a need to investigate the evolution of UHI of mid-size cities (under 1 million inhabitants) like Sari to understand how human activity affects and, more importantly, how the increase in vegetation can play a role in mitigating the urban heat island effects.

Conclusions
The proposed spatiotemporal analysis method, amalgamating Emerging Hot Spot Analysis (EHSA) and the Mann-Kendall test, offers a theoretical framework to enhance our comprehension of the intricate interplay among urbanization, climate change, and the environment.By employing EHSA, the approach identifies spatial clusters with significant changes in vegetation over time, shedding light on the dynamic patterns of urban greenery.Simultaneously, the Mann-Kendall test explores temporal trends, revealing shifts in both vegetation indices and land surface temperatures.This combined methodology allows for a nuanced understanding of how urbanization processes, reflected in changing vegetation, interact with the broader context of climate change.The analysis contributes theoretically by unraveling the spatiotemporal nuances of urban heat island effects, offering valuable insights into the evolving environmental dynamics within urban areas.The interconnectedness of urbanization, climate change, and the environment is thus elucidated through a comprehensive examination of both spatial and temporal dimensions.This study aimed to characterize the urban heat island effects in a mid-size Iranian city representing a humid subtropical climate.Thermal bands and infrared bands of Landsat Areas with higher density (station 2-Local Climate Zone 2) record high UHII day and night (4.0 • C and 4.2 • C).Areas with lower density (station 3-LCZ 5) show high values of UHI during daytime (up to 4.8 • C) and lower values of UHI intensity during nighttime (up to 2.8 • C).

Figure 3
also displays the five-year average dispersion of average LST and UHII.The spread pattern of LST in urban areas is changing in this figure; the city level's LST range was between 29.38 °C and 30.1 °C between 2001 and 2005.The highest LSTs were found in the city center, and the lowest LSTs were in the east and near the Tajen River.The LST distribution pattern in 2006-2010 indicates a range of changes between 29.38 °C and 30.0 °C, like the prior decade.On the other hand, the distribution of the highest LST has been extended to parts of the city's south.The variation in LST increased from 2011 to 2015, reaching 29.58 °C to 30.40 °C.The LST intensity distribution, on the other hand, was like the previous period.The range of LST variations was like the previous period from 2016 to 2020.During this time, the south of the city had the highest LSTs.The LST was at its lowest in 2001-2005 and 2006-2010, with the city's major areas ranging from 29.6-29.51

Figure 2 .
Figure 2. Changes in LST ( • C) (a) and average UHI (b) in urban and rural areas.

Figure 3
Figure 3 also displays the five-year average dispersion of average LST and UHII.The spread pattern of LST in urban areas is changing in this figure; the city level's LST range was between 29.38 • C and 30.1 • C between 2001 and 2005.The highest LSTs were found in the city center, and the lowest LSTs were in the east and near the Tajen River.The LST distribution pattern in 2006-2010 indicates a range of changes between 29.38 • C and 30.0 • C, like the prior decade.On the other hand, the distribution of the highest LST has been extended to parts of the city's south.The variation in LST increased from 2011 to 2015, reaching 29.58 • C to 30.40 • C. The LST intensity distribution, on the other hand, was like the previous period.The range of LST variations was like the previous period from 2016 to 2020.During this time, the south of the city had the highest LSTs.The LST was at its lowest in 2001-2005 and 2006-2010, with the city's major areas ranging from 29.6-29.51• C and 29.81-29.9• C.This range is equivalent to 80.15% of the city's area in the first period and 94.26% of Sari's area in the second period.However, the LST class range increased between 2010 and 2015, and the greatest urban area was between LST classes 29.9-29.81• C and 30.30-30.21• C. So, from 2010 to 2015, the area of coverage of these two tiers covered 97.19% of the city, and from 2016 to 2020, it covered 96.44% of the city.According to Figure 3, the range of LST changes in the city is generally increasing.

Figure 3 .
Figure 3. Area changes related to LST classes based on 5-year averages.Furthermore, Figure 4′s assessment of UHII changes shows that the largest UHII was seen between 2001 and 2005.As a result, the distribution of spots with intensity classes ranging from 0.40-0.31°C to 0.50-0.41°C covered roughly 29.14% of the entire city in the center region, while the scope of these two classes was 16.45% from 2006 to 2010, 20.05% from 2015 to 2011, and 13.13% of the city's total area from 2016 to 2020.The figure shows that the severity of Sari's UHI has worsened over the last two decades.Figure5shows the area changes of the UHII in different classes and the skewness of the UHII trends to be positive in 2001-2005 compared to 2016-2020.In other words, the greater part of city areas with decreased UHII occurred between 2016 and 2020.An important sign of the UHII is that, in addition to existing in areas with low construction density, the UHI has not spread.Rather, the Tajen River's role in weakening the power of the UHI is pretty evident, and it has remained so across all four periods.In general, the results of this phase of the study indicate a decrease in UHII with the rise in NDVI in recent years.These trends are particularly noticeable in the second through fourth decades.The overall average NDVI for these decades is 0.211, 0.228, and 0.238, respectively.Correspondingly, the average UHII for the second and third decades is 0.197 °C, while for the fourth decade, it is 0.172 °C.Conversely, the findings reveal an increasing trend in both high vegetation and low vegetation areas from the second to the fourth period.This growth in green spaces within the urban environment directly contributes to the reduction in UHII.

3 .
Area changes related to LST classes based on 5-year averages.Furthermore, Figure 4's assessment of UHII changes shows that the largest UHII was seen between 2001 and 2005.As a result, the distribution of spots with intensity classes ranging from 0.40-0.31• C to 0.50-0.41• C covered roughly 29.14% of the entire city in the center region, while the scope of these two classes was 16.45% from 2006 to 2010, 20.05% from 2015 to 2011, and 13.13% of the city's total area from 2016 to 2020.The figure shows that the severity of Sari's UHI has worsened over the last two decades.Figure 5 shows the area changes of the UHII in different classes and the skewness of the UHII trends to be positive in 2001-2005 compared to 2016-2020.

Figure 4 .
Figure 4. Spatial-temporal changes of LST, UHI and NDVI based on 5-year average intervals.Figure 4. Spatial-temporal changes of LST, UHI and NDVI based on 5-year average intervals.

Figure 4 .
Figure 4. Spatial-temporal changes of LST, UHI and NDVI based on 5-year average intervals.Figure 4. Spatial-temporal changes of LST, UHI and NDVI based on 5-year average intervals.

Figure 5 .
Figure 5. Area changes related to UHII classes based on 5-year averages.

Figure 6 .
Figure 6.Changing trend map of Sari's UHII calculated using the Mann-Kendall method .

Figure 5 .
Figure 5. Area changes related to UHII classes based on 5-year averages.

Figure 5 .
Figure 5. Area changes related to UHII classes based on 5-year averages.

Figure 6 .
Figure 6.Changing trend map of Sari's UHII calculated using the Mann-Kendall method .

Figure 6 .
Figure 6.Changing trend map of Sari's UHII calculated using the Mann-Kendall method.

Figure 7 .
Figure 7. Changing trend map of Sari's NDVI calculated using the Mann-Kendall method .

Figure 7 .
Figure 7. Changing trend map of Sari's NDVI calculated using the Mann-Kendall method.

Figure 8 .
Figure 8. Area changes related to NDVI classes based on 5-year averages.

Figure 8 .
Figure 8. Area changes related to NDVI classes based on 5-year averages.

Figure 9 .
Figure 9. Revealing spatiotemporal anomalies of UHII in Sari City using EHSA analysis.

Figure 9 .
Figure 9. Revealing spatiotemporal anomalies of UHII in Sari City using EHSA analysis.
7 and 8 satellites were used to generate a dataset between 2000 and 2022 for Sari, Iran.The Mann-Kendall and Emerging Hot Spot Analysis techniques were used as analytical methods for spatiotemporal heat island effect analysis.The results indicate that the LST of Sari has been increasing in recent years due to global warming.Other variables, such as human indicators (AI), could have influenced it.The urban LST area variations from 2011 to 2015 were found in the largest part of Sari City in the LST class of 30.2-30.11• C, which equates to 35.11%.It is like 39.46% of the overall city in the LST class of 29.8-29.71• C from 2006 to 2010.The Normalized Difference Vegetation Index increase over the study period decreased the UHI intensity.The city's LST increased by approximately 0.30 • C between 2011 and 2022.However, the intensity of the UHII decreased during the study period.An average −0.031 • C/decade decrease in the UHII has been experienced.The findings prove the positive role of urban green vegetation and its influence in reducing heat stress.Future research should investigate further effective measures to introduce nature-based solutions concerning water and energy use in urban dense agglomerations.Also, field surveys and local measurement campaigns are needed to characterize the spatial and temporal variability of the UHI.

Table 1 .
The frequency of LST days was calculated from Landsat 7 and 8 satellite images.

Table 1 .
The frequency of LST days was calculated from Landsat 7 and 8 satellite images.

Table 2 .
Urban vegetation classes and NDVI value.

Table 4 .
Spatial distribution of the significant level of change in the UHII of Sari City by Mann-Kendall method.

Table 5 .
Spatial distribution of the significant level of change in the NDVI of Sari City by Mann-Kendall method.

Table 6 .
The area related to UHII anomalies of Sari City according to each of the EHSA patterns.

Table 6 .
The area related to UHII anomalies of Sari City according to each of the EHSA patterns.
• C compared to 2001 to 2010, according to an analysis of the Sari urban LST map.This amount is significant for 20 years in a town of moderate size like Sari.Based on an investigation of urban LST area variations from 2011 to 2015, the largest part of Sari City in the LST class of 30.2-30.11• C equates to 35.11%.At the same time, it is 39.46% of the overall city in the LST class of 29.8-29.71• C from 2006 to 2010.