Synergies between Urban Heat Island and Urban Heat Wave Effects in 9 Global Mega-Regions from 2003 to 2020

Abstract

Global urbanization significantly impacts the thermal environment in urban areas, yet urban heat island (UHI) and urban heat wave (UHW) studies at the mega-region scale have been rare, and the impact study of urbanization is still lacking. In this study, the MODIS land surface temperature (LST) product was used to depict the UHI and UHW in nine mega-regions globally between 2003 and 2020. The absolute and percentile-based UHW thresholds were adopted for both daily and three-day windows to analyze heat wave frequency, and UHW magnitude as well as frequency were compared with UHI variability. Results showed that a 10% increase in urban built-up density led to a 0.20 °C to 0.95 °C increase in LST, a 0.59% to 7.17% increase in hot day frequency, as well as a 0.08% to 0.95% increase in heat wave number. Meanwhile, a 1 °C increase in UHI intensity (the LST differences between the built-up and Non-built-up areas) led to a 2.04% to 92.15% increase in hot day frequency, where daytime LST exceeds 35 °C and nighttime LST exceeds 25 °C, as well as a 3.30% to 33.67% increase in heat wave number, which is defined as at least three consecutive days when daily maximum temperature exceeds the climatological threshold. In addition, the increasing rates of UHW magnitudes were much faster than the expansion rates of built-up areas. In the mega-regions of Boston, Tokyo, São Paulo, and Mexico City in particular, the increasing rates of UHW hotspot magnitudes were over 2 times larger than those of built-up areas. This indicated that the high temperature extremes, represented by the increase in UHW frequency and magnitudes, were concurrent with an increase in UHI under the context of climate change. This study may be beneficial for future research of the underlying physical mechanisms on urban heat environment at the mega-region scale.

1. Introduction

It is projected that the population in the urban areas will reach 5 billion by 2030 [1]. Accelerating global urbanization processes have led to the agglomeration of adjacent urban areas known as mega-regions. Numerous land surface modification has also triggered the conversion of vegetated surfaces to buildings or impervious surfaces [2,3,4]. This has exerted increasingly stronger impacts on the ecosystem, biodiversity, and local climate, and has brought about many negative environmental repercussions, such as urban heat phenomena [5]. Specifically, urban heat island (UHI) phenomena, where urban areas tend to have a higher near-surface air temperature (NSAT) or land surface temperature (LST) relative to the surrounding hinterland [6,7,8], have been widely used to demonstrate the magnitude of urban heat phenomena. Meanwhile, urban heat wave (UHW) phenomena, defined as events associated with particularly high and sustained NSAT, have also been used to quantify the frequency, intensity, and duration of urban heat events [9,10,11,12]. These urban heat phenomena have a large impact on human life, health, and well-being in a variety of ways. Consequently, the understanding of urbanization processes and the resulting urban heat environment—both their patterns and intensity—are of particular importance to natural resource use, public health, and global environmental changes [13,14,15].

NSAT is normally measured at thermometer shelters 1.5–2 m above the ground at weather stations with high accuracy and temporal resolution [16,17]. However, these data are heavily dependent on the regional infrastructure for weather data collection. They are collected as point samples whose distribution is rarely designed to capture the range of climate variability within a region. The spatial information available on air temperature is thus often limited, especially in areas with rugged terrain and harsh environments. The need for spatial information and real-time access to data has thus driven many researchers to look for satellite-based methods. Meanwhile, LST is instead retrieved from thermal infrared or microwave spectral measurements. It is available from spaceborne radiometers and provides a synoptic spatial perspective that can be very important in the modeling and monitoring of UHI processes at different spatial scales [18]. With advances in technology, communications, and computing power, data from satellite sensors have become easily accessible. In particular, with the advantages of wide coverage and short revisit intervals, various thermal infrared remote sensors, such as the Advanced Very High Resolution Radiometer (AVHRR) [19,20], Moderate Resolution Imaging Spectroradiometer (MODIS) [21,22], Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [23,24], and Landsat Thematic Mapper [25,26], have been used to measure LST intensity to quantify urban heat phenomena. Generally, though the NSAT and LST feature different measurement techniques, diurnal phases, and responses to atmospheric conditions [27], LST is strongly consistent with NSAT. Many previous studies have also retrieved high-solution NSAT data from MODIS products over different ecosystems, and the difference between LST-retrieved NSAT and LST ranges from 0~1.73 °C [28,29], implying that LST can also provide a reliable measure of the near-surface thermal state that is spatially distributed and temporally consistent [30,31].

However, most of the satellite-based urban heat phenomenon studies have focused on the magnitude and intensity of UHI based on LST gradient [32,33,34], while LST data have rarely been used to quantify the frequency and duration of urban heat events. One reason is that measuring the duration or frequency of urban heat events normally requires consecutive availability of NSAT data at a sub-daily scale to depict diurnal temperature dynamics. As a result, the utility of 90–100 m Landsat/ASTER thermal data with 16-day temporal resolution is limited in detecting urban heat frequency or intensity without a series of downscaling [35], sharpening [36], image fusion [37], and disaggregation [38] methods to enhance the temporal resolution. By contrast, MODIS LST data provide great potential with 1 km spatial resolution and 12 h temporal resolution, but also pose considerable challenges in computational burden. The challenge may be met with the assistance of cloud-based computation platforms such as the Google Earth Engine (GEE) [39], which has a large volume of global satellite imagery storage and high computing capability. This may dramatically improve the quantitative analysis and understanding of the frequency and intensity of urban thermal phenomena.

Therefore, the first research goal in this study is to apply the high spatial resolution and moderate temporal resolution MODIS LST data to measure the frequency and intensity of urban thermal events based on the GEE platform. As there is a strong correlation between LST and NSAT, we used LST as a surrogate of NSAT to depict UHW using a similar definition. To distinguish it from traditional NSAT-based UHW delineation, the LST-based UHW evaluation is referred to as the pixel-based UHW evaluation in this study. The definitions of UHW in previous studies include a number of percentile-based, absolute, duration-based, or range-based NSAT thresholds [10,40,41,42,43], which are largely subject to the context and location. For example, aside from the 17 extreme temperature indices listed in Climate Change Detection and Indices (ETCCDI) [44], more than 20 other UHI/UHW indices exist, which were specifically designed to depict certain impact communities (e.g., human mortality, wildlife, agriculture, and transport). The robustness and transferability of these indices have not been verified in many areas due to the lack of high-quality observations [45]. Different indices were adopted in these studies, making the inter-comparison between urban heat phenomena studies difficult [46]. Therefore, there is an urgent need for integrating remote sensing data sets based on an empirically verified theoretical foundation upon which general indices of LST-based UHI/UHW are applicable and transferable to different impact communities or study areas. The theoretical framework will potentially enable comparisons of past, current, and future urban heat phenomena, as well as the spatial magnitude, temporal intensity, and frequency of different LST-based urban heat scenarios.

At the same time, previous satellite-based analyses indicated that different land use and land cover characteristics, such as impervious surface area (ISA), urban green space coverage, and landscape diversity of urban metropolis [47,48,49,50,51], are the most influential factors on the magnitude of UHI [52,53,54]. In particular, ISA, defined as anthropogenic features through which water cannot infiltrate into the soil (e.g., roads, driveways, side-walks, parking lots, rooftops, etc.), has proved to be the primary driving factor for the increase in temperature [55,56,57]. Previous studies have also shown that ISA accounted for over 60% of the total LST variance of urban settlements within forests at mid-to-high latitudes globally [58,59]. The relationship between LST and ISA has been used as a potentially powerful tool for quantifying the contribution of land cover changes induced by human activities to urban heat phenomena [60,61,62].

As the ongoing global urban expansion strongly directs the evolution of city networks into a new spatial scale—the mega-region scale [30,31]—urban heat phenomena represented by LST features have been shown to be more pronounced in the ISA. For example, the generally quasi-linear relationship between LST and the built-up density (depicted by ISA fraction) has been examined in the most densely populated cities in the United States [62] and China [63], but very few studies have focused on the spatiotemporal trends of urban heat intensity and frequency, as well as their relationship with ISA fractions at the mega-region scale. It is worth exploring whether the relationship between urban heat phenomena and built-up density adopts a similar quasi-linear or a non-linear pattern. Therefore, there is a strong impetus to systematically evaluate the spatiotemporal trends of urban heat phenomena with ISA fractions in different mega-regions [64,65,66,67].

This study aims to analyze the spatiotemporal patterns of urban heat intensity and frequency based on LST data from MODIS in 9 mega-regions globally between 2003 and 2020. Further, the influence of ISA fractions on UHI/UHW is evaluated using an urban built-up density product [68]. Two key research questions addressed are:

1.

What is the added value of coupling UHI and UHW analysis in characterizing urban heat phenomena using remote sensing data?

2.

What is the influence of ISA fractions on UHI/UHW phenomena?

This study is among the first studies to evaluate urban heat frequency in depth with remotely sensed LST data using indicators originally developed for NSAT. We specifically focus on exploring the spatiotemporal interaction of urban heat intensity and frequency patterns, as well as its relationship with built-up density in representative mega-regions globally.

2. Study Area and Data

2.1. Study Area

The magnitude of population (over 10 million) and geographical expansion are two typical defining features of mega-regions [69,70,71], but the delineation of the geographical boundary of a mega-region is still vague. Based on a comprehensive literature review, nine representative mega-regions across the globe with a population larger than 20 million by 2000 were selected. Out of the nine mega-regions, three are located in coastal deltas; we thus named them after their respective deltas, the Yangtze River Delta (YRD), the Pearl River Delta (PRD), and the Nile Delta. The other six mega-regions contain mega-cities with a population over 10 million; we therefore used the mega-city‘s name to refer to the mega-region, including Paris, Mexico City, Tokyo (mainly from Tokyo to Osaka), São Paulo (mainly from São Paulo to Rio de Janeiro), Jakarta (mainly from Jakarta to Surabaya), and Boston (mainly from Boston to Washington). The province-scale geospatial data from the GADM dataset (https://gadm.org/ (access on 5 July 2021)) were used to determine the geographical boundaries of the nine mega-regions (Figure 1). These regions feature different climate types, including tropical climate (Nile Delta, Jakarta), subtropical climate (PRD, São Paulo, Mexico City), and temperate climate (YRD, Boston, Tokyo, Paris), with the mean summer temperature (estimated from MODIS LST product (MYD11A1)) ranging from 19.31 °C to 33.17 °C. Based on the urban built-up area data provided by Gong et al. [68], the highest growth rate of built-up density between 2000 and 2018 was seen in the mega-region of YRD (53.56%), followed by the mega-regions in the Nile Delta (31.45%), Paris (28.43%), PRD (18.51%), Boston (14.87%), Mexico City (13.44%), Jakarta (13.15%), São Paulo (5.36%), and Tokyo (4.02%). The highest growth rate of population density between 2000 and 2015 was in the mega-region of PRD (44.50%), followed by the Nile Delta (31.76%), Mexico City (25.76%), YRD (23.88%), Jakarta (20.82%), São Paulo (18.41%), Paris (10.36%), Boston (8.44%), and Tokyo (0.96%), based on the geostatistical data from Gridded Population of the World (V4) [72]. More details are shown in Table S1 in the Supplementary Materials.

2.2. Data

In this study, the LST product (MYD11A1) from MODIS Aqua [73,74] between 2003 and 2020 under clear-sky (99% confidence) condition was used. The data were provided daily at 1 km resolution. The quality control procedures included removing cloudy images using the quality control band of MOD11A1 data to only keep the pixel with quality flag = 0 (good data quality) and average emissivity flag = 0 (average emissivity error ≤ 0.01). The LST observations at the overpass time of 13:30 and 01:30 local time were taken as proxies of the daytime and nighttime temperature. The data accuracy was shown to be better than 1 K for most cases, and the difference with in situ measurements was generally less than 5% in urban areas [75,76]. Pixels with water bodies, or those with no value or of bad quality, were masked out through quality control prior to evaluation. The seasonal UHI and UHW features were then calculated for the summer months between 2003 and 2020 for each mega-region, respectively. Here, summer months were defined as December to February of the following year in the mega-region of São Paulo, and as June to August in the other mega-regions.

Urban built-up area data were provided by Gong et al. [68] with annual maps of global artificial impervious areas between 2003 to 2018. This product is derived at 30 m resolution from cloud-free Landsat TM images, and the classification accuracy is generally higher than 0.9 for most mega-regions. Here, the data of 2018 were assumed representative for 2019 and 2020. More details about the built-up area product can be found in Gong et al. [68].

3. Methodology

This study measures the urban heat environment based on LST data from the MODIS daily composite product in nine mega-regions between 2003 and 2020. The urban heat environment was evaluated using four indicators (

Table 1. Indicators used to characterize UHI and UHW patterns.

	Feature	Indicator	Definition	Reference
UHI	Intensity	mean LST	mean LST in the urban areas of each mega-region	[26,77]
		LST standard deviation	LST standard deviation in the urban areas of each mega-region
UHW	Frequency	Combined hot days and tropical nights (CHT)	Seasonal occurrence of summer days with maximum temperatures exceeding 35 °C and minimum temperatures exceeding 20 °C	[10,31,41,43,78]
		Heat wave Number (HWN)	A heat wave event is defined as at least 3 consecutive days when daily maximum temperature exceeds the climatological 90th/97.5th percentile, and is anomalously warm compared to the historical mean LST.

). These indicators were designed to measure the intensity and frequency of UHI and UHW, and aimed (1) to evaluate a number of pixel-based urban heat indicators on their feasibility across different climates and (2) to facilitate a spatially consistent comparison of patterns and trends of observed thermal phenomena among different mega-regions.

3.1. Urban Heat Island Evaluation

In order to demonstrate the UHI and UHW trends induced by rapid urbanization, we first aggregated the 30 m resolution built-up area data to 1 km resolution and calculated an urban built-up density index (UDI), which was then used to delineate built-up areas (UDI > 0). The magnitude of UHI was defined as the total built-up areas, and the intensity of UHI was evaluated by comparing the LST differences in the built-up areas and Non-built-up areas (UDI = 0). The UDI data were further categorized into 10 levels, with UDI values ranging from 0 to 1 with an increment of 0.1. The mean and standard deviation (std) of LST, as well as the LST correlation with UDI, were calculated for each UDI level. The formulas of UDI and UHI are shown as follows:

$${\mathit{UDI}}_{\mathit{pixel-based}} = \sum IS{A}_{\left(30m\right)}/\sum Are{a}_{\left(1km\right)}$$

(1)

UHI = T_{(UDI > 0)} − T_{(UDI = 0)}

(2)

3.2. Urban Heat Wave Evaluation

The analysis of UHW intensity and magnitude depends to a large extent on the definition of heat waves. Here, two indicators were adopted to characterize heat wave events: the combined hot days and tropical nights (CHT), and the heat wave number (HWN). The CHT indicator is based on an absolute LST threshold, which represents the average number of days with daytime temperature exceeding 35 °C and nighttime temperature exceeding 20 °C [12]. The HWN, on the other hand, is a relative indicator that measures heat wave intensity, and depends on the temperature threshold determined by the LST climatology of the study area. Here, a heat wave event is defined as a spell of at least three consecutive days with daytime temperature exceeding the local 90th percentile or 97.5th percentile of the study period [78]. The HWN thresholds were determined by:

T_97.5 = T_{seasonal_mean} + 2 T_{seasonal_std}

(3)

T₉₀ = T_{seasonal_mean} + 1.28 · T_{seasonal_std}

(4)

where Tseasonal_mean and Tseasonal_std are the historical (2003–2020) average LST and standard deviation of each pixel. T₉₀ and T_97.5 are the 90th and 97.5th percentiles of the LST climatology of each mega-region from 2003 to 2020. To eliminate the impact of cloud contamination on the LST data availability, the derived CHT and HWN were normalized by the number of valid observations per season for each pixel. The UHW intensity was defined as the value of CHT and HWN, and the UHW magnitude was defined as the area with at least one heat wave event.

3.3. The Hotspot of Urban Thermal Environment

Percolation theory was originally developed in statistical physics and mathematics to model the emergent structures of urban clusters on a random graph [79]. Urban researchers have used percolation theory to model urban boundaries and to understand the critical phenomena of cities [78,80]. We applied percolation theory to delimit hotspots of urban thermal environment using LST, CHT, HWN90th, and HWN97.5th, separately. As shown in Figure 2, for an L × L lattice, each pixel can be occupied with probability p (e.g., LST, CHT, or HWN), and adjacent occupied pixels form a cluster. The size of the largest cluster remains stable for small p values. As p becomes larger, the size of the largest cluster starts to increase rapidly and gradually becomes stable again at a certain point, following a power-law distribution. This point is called the critical point, or the percolation threshold. Here, we consider the threshold at this critical point to be the optimal threshold. Once the optimal threshold was determined, we were able to delineate the hotspots of urban thermal environment based on different indicators.

4. Results

4.1. The Trends of Urban Heat Phenomena in the Built-Up and Non-Built-Up Areas

The temporal dynamics and variability of the mean LST in the study area from 2003 to 2020 were first assessed (Figure 3). Linear regression showed an increasing trend in the mean LST in eight mega-regions except Jakarta. An increasing trend was also found in the UHI intensity (mean LST difference between the built-up and Non-built-up areas) in all nine mega-regions, which was significant (p < 0.05) in all regions except for Mexico City and Paris. The largest annual average UHI intensity was seen in Tokyo (6.59 °C), followed by the mega-regions of PRD (3.64 °C), São Paulo (3.23 °C), Boston (3.19 °C), Jakarta (3.12 °C), YRD (2.70 °C), Mexico City (1.94 °C), and Paris (1.37 °C). Meanwhile, the largest increasing rate of UHI intensity was seen in Paris (11.35%), followed by YRD (10.84%), Nile (9.10%), PRD (8.83%), Jakarta (5.29%), Tokyo (4.04%), São Paulo (1.63%), Boston (1.15%), and Mexico City (0.79%). The UHI intensity in the summer daytime was generally higher than that in the summer nighttime in these mega-regions (shown in Figure S1 in the supplementary material). This result is consistent with previous research, such as Ayansina [81] and Zhou [82]. The UHI intensity can be attributed to the varying mechanisms underlying the daytime and nighttime UHI effects. Within the mega-region, the densely vegetated areas which are mainly located in the area with (UDI = 0) play a crucial role in cooling the earth surface via transpiration during the daytime. Comparatively, vegetation demonstrates a weak effect over nighttime because transpiration is significantly reduced at night [82]. Consequently, the UHI intensity (i.e., LST difference between impervious surface areas and non-impervious surface areas) during the daytime are normally larger than that during the nighttime.

Different from the other mega-regions, the Nile mega-region exhibited an evident ‘cold island effect’ (UHI intensity < 0), where the urban area has a lower temperature than the surrounding areas. Previous studies have also reported the ‘cold island effect’ in some arid and semi-arid environments, such as in Athens [83,84]. This may be caused by the fact that the major land cover types in this mega-region are bare soil and sparse vegetation [85,86]. Based on the MODIS yearly land use product (MCD12Q1.006), over 75% of the area is covered by bare soil from 2003 to 2020. During the daytime, the open plain mainly covered by sparse vegetation or bare soil becomes heated faster than the densely built-up areas, which are mainly covered by materials with higher thermal inertia, such as concrete and asphalt [87]. In addition, bare soil has a lower albedo (typical value around 0.17) than concrete (typical value around 0.55), leading to higher net radiation in the Non-built-up areas. A dense building envelope with specific orientation, or building external walls with high-reflectivity materials in hot-arid climates, can also reduce heat gains from solar radiation and create a cooler environment through reducing cooling load accordingly [88]. During the nighttime, urban areas tend to be cooler as well, as a result of the lower thermal inertia of the soil compared to concrete [8,89].

The mega-region of Mexico City also demonstrated a ‘cold island effect’ in the night-time. The region features a distinct wet season and dry season, and temperatures remain fairly constant throughout the year. The wet season lasts from May to October and features heavy rainfalls. The ‘cold island effect’ in the nighttime is possibly induced by the relatively high vegetation cover and soil moisture in the wetland and nearby buffer zones in the rainy season. Similar nighttime heat island phenomena were also found in Mexico City in previous research [90].

The mean CHT trends in the nine mega-regions are shown in Figure 3. Out of the nine mega- regions, eight (except for Paris) showed a positive CHT trend, and the trends were significant in three mega-regions (Nile, YRD and São Paulo). The largest increase in CHT intensity between 2003 and 2020 was seen in the Nile mega-region (126.27%), followed by the mega-regions in YRD (88.11%), São Paulo (73.46%), Tokyo (49.73%), Jakarta (19.86%), PRD (12.97%), Boston (4.12%), and Mexico City (3.52%). It is worth pointing out that the CHT value in 2003 in Paris was significantly higher than that from other years, as a result of an extreme heat wave event that occurred in Western Europe in 2003. If the 2003 value is excluded, the CHT intensity in Paris would also demonstrate a significant positive trend, with an increase of 11.95% from 2004 to 2020 in the built-up area. In addition, except for the mega-region of Mexico City, the increasing rate of CHT intensity in the built-up area was also higher than that in the Non-built-up areas.

The trends of HWN90th are shown in the third column of Figure 3. Similarly to the trends of CHT, five out of the nine mega-regions (Nile, Tokyo, Paris, YRD, and São Paulo) also demonstrated an increasing trend of HWN90th intensity from 2003 to 2020 in the built-up areas, but only one of them was significant (Nile). The largest increase in HWN90th intensity was seen in the Nile mega-region (13.24%), followed by the mega-regions in YRD (7.58%), Tokyo (5.01%), São Paulo (1.92%), and Paris (0.01%). The HWN90th increased more rapidly in the built-up areas than in the Non-built-up areas. Though HWN90th fluctuated in the mega-regions of Boston, PRD, Jakarta, and Mexico City over the study period, the HWN90th values in the built-up areas were almost always higher than those in the Non-built-up areas. Nile is the only mega-region with higher HWN90th values in the Non-built-up areas, which is consistent with the trend of the UHI intensity. The temporal trends of HWN97.5th are consistent with the results using HWN90th (shown in the fourth column of Figure 3).

4.2. Influence of Built-Up Density on UHI/UHW Phenomena

To further evaluate the influence of built-up density on the urban heat phenomena, the mean LST within the mega-regions was calculated at 10 levels with UDI from 0 to 1 with an increment of 0.1 (e.g., Level 1 indicates areas with 0 < UDI < 0.1). The box plot of mean LST within each UDI level during the study period is shown in Figure 4. The Pearson correlation coefficient (R²) was used to quantify the strength of the relationship between UDI and LST. A positive correlation between UDI and the mean LST was observed in all mega-regions. The R² ranged between 0.1 and 0.8 and was significant (p < 0.01) in all mega-regions except for Nile. It is worth noting that the UDI ~ LST relationship was still positive in the urban areas in the mega-region of Nile, despite the ‘cold island effect’ observed.

The mean CHT and the mean HWN90th intensity for different UDI levels are shown in the second and the third column of Figure 4. A consistent positive correlation (R² = 0.47~0.99, p < 0.01) between UDI and CHT is seen in all the mega-regions from 2003 to 2020, except for the mega-region of Mexico City. Moreover, within these positive correlations, 118 out of the 162 relationships between UDI and CHT were significant. The HWN90th metric also exhibited a significant positive (R² > 0.72, p < 0.01) relationship with UDI in all mega-regions except for Nile and Mexico City. This may relate to the influence of bare soil and wetland in these two regions. In addition, the correlations between UDI and HWN90th are generally higher than those between UDI and HWN97.5th.

Overall, the increase of 10% in UDI corresponds to around a 0.20 °C to 0.95 °C increase in LST, a 0.59% to 7.17% increase in CHT, as well as a 0.08% to 0.95% increase in HWN90th. The positive correlations between UDI and LST/CHT/HWN indicate that the built-up areas realistically represent the compound UHI and UHW effects. In addition, when taking the four types of climate zones among the mega-regions: (1) Humid subtropical (YRD, PRD, Tokyo, Boston, Mexico City, São Paulo); (2) Tropical rainforest (Jakarta); (3) Marine west coast climate (Paris); (4) Desert climate (Nile), into consideration, we observed some patterns potentially associated with climatic conditions. For example, the mega-regions in the humid subtropics more often exhibit stronger UHI intensity (LST difference between built-up and Non-built-up areas) [91,92]. In contrast, thanks to the moderating effect of the ocean [93], the difference between thermal environments of built-up and Non-built-up areas is less evident. The desert climate of the Nile delta on the other hand, contributes to the unique ‘cold island effect’, which is the result of mixed influence of both the climate and the landscape [94].

4.3. The Spatiotemporal Patterns of Urban Heat Phenomena

The mean of LST, CHT, and HWN intensities during the study period are shown in Figure 5. The LST, CHT, and HWN spatial patterns were generally in line with previous studies based on in situ data or model simulations [95]. Areas with high CHT and HWN intensities were predominantly located in the urban areas with high built-up density. The core urban areas have marked influence on the aggregation of urban heat wave effects, and the larger the core areas, the more significant the urban heat aggregation effect. In particular, some urban agglomerations in the mega-regions of Boston, Tokyo, Jakarta, YRD, São Paulo, and Mexico City have merged and evolved into urban corridors along the coastal areas, where the CHT and HWN intensities demonstrated a declining trend from the multiple cores of these urban corridors to the fringe areas. In addition, though the core urban areas of Paris were not located in the coastal area, its urban heat phenomenon pattern also revealed a declining trend from the city center to the fringe areas.

When the increasing rates of CHT magnitudes (the proportion of built-up areas with at least one CHT event) were taken into consideration, the increasing rates of CHT magnitudes in six mega-regions were much higher than the expansion rate of built-up areas. In the mega-regions of Boston, Tokyo, São Paulo, and Mexico City in particular, the increasing rates of CHT magnitude were more than 2 times higher than those of built-up areas. However, the higher increasing rates of UHW phenomena than those of built-up areas were not observed using the HWN indicators, except for in Paris.

The percolation-based hotspot detection result is shown in Figure 6. Different indicators demonstrated different UHI and UHW clustering patterns in the nine mega-regions. As shown in

Table 2. The largest proportion of UHI hotspot magnitude based on different UHI and UHW patterns.

	LST	CHT	HWN90th
YRD	97.32%	30.59%	6.14%
Boston	97.22%	18.40%	61.34%
PRD	95.97%	29.08%	1.18%
Nile	95.67%	86.34%	80.64%
Mexico City	93.25%	10.05%	51.83%
Paris	93.06%	4.93%	27.80%
São Paulo	91.69%	46.05%	75.94%
Tokyo	88.13%	2.22%	11.89%
Jakarta	83.82%	17.76%	15.59%

, the largest proportion of UHI hotspot magnitude (the area of the largest cluster/the area of the mega-region) based on LST was YRD, followed by Boston, PRD, Nile, Mexico City, Paris, São Paulo, Tokyo, and Jakarta, while the largest proportion of CHT-based hotspot magnitude was Nile, followed by São Paulo, YRD, PRD, Boston, Jakarta, Mexico City, Paris, and Tokyo. In terms of HWN90th, the largest hotspot magnitude proportion was also found in Nile, followed by São Paulo, Boston, Mexico City, Paris, Jakarta, Tokyo, YRD, and PRD. The spatial extent of the UHI and UHW clusters detected using different indicators was mostly consistent in the mega-regions of Boston, São Paulo, and Nile. In other mega-regions, such as Jakarta, Tokyo, Paris, PRD, and YRD, the clusters aggregated more evidently over the built-up areas, despite the dissimilar hotspot patterns using different indicators.

The study period was further divided into three 6-year periods to analyze the percolation-based hotspot patterns (more details are included in Tables S6–S9 in the Supplementary Materials). We found that the average LST window width (as illustrated in Figure 2) steadily increased in Mexico City (by 3 °C), Paris (by 2 °C), and PRD (by 2 °C), while it fluctuated in the other mega-regions. Similarly, the CHT-based hotspots also demonstrated a steady increase in six mega-regions, with the largest increase occurring in Boston (28%) and the lowest occurring in Mexico City (7%). Though the HWN-based hotspots did not demonstrate significant increasing trends within the three periods, the overall patterns of UHI/CHT/HWN-based hotspots reveal that anthropogenic activities have detectable influence on the thermal environment at the mega-region scale and should be considered in future research which explores the mechanistic drivers of heat waves.

5. Discussion

5.1. The Coupling Effect of UHW and UHI Phenomena

The correlation analysis between daytime/night-time mean LST, CHT, HWN90th and HWN97.5th, and the global mean temperature (gridded data from 2003 to 2018 provided by NOAA, https://www.ncdc.noaa.gov/temp-and-precip/global-maps/ (accessed on 10 March 2020)) was performed (shown in Table S10 in the Supplementary Materials). Except for the mega-regions, it is shown that the global mean temperature during the summer season increased by 0.34 °C from 2003 to 2018. The first peak was in 2005, with the temperature having grown by 0.15 °C from 2003. The second peak was in 2010 with the temperature having grown by 0.19 °C from 2005, and the third peak was in 2016 with the temperature having grown by 0.15 °C from 2010. The increase in global mean temperature coincides with the occurrence of natural phenomena such as ENSO or volcanism, which are known to affect the duration of heat waves [96]. Following the increase in global mean temperature, the frequency and intensity of the temperature extremes evaluated by multiple indicators (including mean LST, HWN97.5th, HWN90th) also increased in most of the mega-regions (except for Nile and Mexico City). All these indicators in the built-up areas demonstrated a more significant trend than that in the Non-built-up areas in correlation with the global average temperature (R² > 0.2, p < 0.01).

Many previous studies in climatology have also demonstrated that small changes in the mean temperature can result in disproportionally large changes in the intensity and frequency of extreme heat events. To evaluate the regional differences in the impact of LST variation, the temporal dynamics of mean LST, CHT, HWN90th, and HWN97.5th are plotted in Figure 7. We calculated the pair-wise correlations between the indicators in each mega-region (90 pairs in total). All the correlations were positive, except that of between daytime and nighttime mean LST in Jakarta, and only two correlations were insignificant. The variations of CHT and HWN97.5th agree very well with the mean LST dynamics, which is consistent with earlier findings that the increase in seasonal average temperature over most regions globally can lead to increased occurrence of temperature extremes [97]. Overall, a 1 °C increase in LST corresponded to a 2.04% to 92.15% increase in CHT. The highest CHT variability was found in the mega-region of Mexico City, and the lowest was in the mega-region of PRD. Meanwhile, a 1 °C increase in LST led to a 3.30% to 33.67% increase in HWN90th. The highest HWN90th variability was found in the mega-region of Paris, while the lowest was in the mega-region of PRD. In particular, the CHT and HWN90th variability in the built-up areas was more dramatic than that in the Non-built-up areas, except for the mega-region of Paris. This indicates that when the increased temperature variability is concurrent with an increase in the average temperature, the high temperature extremes represented by the increase in UHI magnitude and UHW frequency would occur more frequently under global warming. Similar findings have also been reported in Diffenbaugh and Giorgi [98] based on CMIP5 data (Coupled Model Inter-comparison Project, Phase 5).

5.2. The Uncertainty of Remote Sensing-Based UHW Measurement

The absolute threshold-based CHT is a transferable indicator that incorporates daily maximum and minimum LST information for UHW delineation, and it provides a spatially uniform measure of heat wave patterns over the nine mega-regions. However, CHT may not be suitable for depicting the spatial-temporal characteristics of heat wave events in regions within the Marine west coast climate zone where the maximum temperature rarely exceeds 35 °C, such as part of the Paris mega-region. For such regions, the percentile-based HWN may be better suited to capture local UHW events, since it takes into account the climatic differences between regions. HWN can also be further analyzed to measure the amplitude, duration, and persistence of UHW events [41,78].

Nevertheless, percentile-based indicators such as HWN are influenced by the time window adopted, as well as the percentile threshold selected. Expert Team on Climate Change Detection and Indices (ETCCDI) have used several indices, such as TN90p, TN95p, TN97.5p, to measure the percentage of days when the daily minimum temperature is above the corresponding percentile [99,100]. These indices have been used in studies based on both observations and model simulations in different regions [41,101]. Meanwhile, previous studies have defined a heat wave event as 3 days [102], or 5 days [100], or 7 days [103], or 15 days [102] that consecutively meet the pre-defined standards. In this study, a heat wave standard similar to the one used in Perkins and Alexander [100] was adopted, and the 90th and 97.5th percentiles of mean LST climatology were selected as the threshold [6,102], respectively. Using the two thresholds, the average number of heat wave events during the daytime ranges from 1.19% to 34.34%, and from 0.04% to 19.78% over the urban areas, respectively. By contract, using a lower threshold (e.g., 80th percentile) or longer time periods (e.g., 6 days/15 days) would result in too many or too few events being detected. The determination of the thresholds is thus a trade-off between effectively detecting representative UHW events and highlighting the severity of extreme UHW events.

Despite the similarities, it is worth noting that the spatial patterns of HWN do not fully coincide with those of CHT. For example, the CHT stands out in the northern part of Mexico City mega-region, which is not seen in HWN patterns. The HWN also does not show the aggregated effect seen in CHT in the mega-region of PRD. This is on one hand caused by the difference of the two indicators by definition, and on the other hand caused by the data availability of the MODIS LST product. The uncertainty of the pixel-based UHI and UHW characteristics in this study is closely related to the amount of available MODIS LST data, which is directly influenced by the weather conditions. The histogram of LST data availability in the built-up and Non-built-up areas is shown in Figure S2 in the supplementary materials. Since the UHW analysis requires consecutive revisits within a certain time period (3 consecutive days in this study), its estimation is not only affected by the magnitude of LST (as in the UHI analysis) but also by the temporal availability of valid LST observations. As a result, the uncertainty of UHW analysis is expected to be larger compared to UHI analysis. In order to reduce the uncertainty, the calculated CHT and HWN values were normalized by the total number of days with valid LST observations to enhance the robustness of the evaluation. Here an implicit assumption is that there are no UHW events under cloudy-sky conditions (i.e., when no LST data are available). Clouds strongly reflect the shortwave radiation from the sunlight, thus greatly reducing the solar radiation that reaches the ground, which essentially limits the magnitude of the surface temperature and subsequently reduces the possibility of heat wave events [104,105].

The earth surface thermal infrared (TIR) radiation conspicuously shows an anisotropic behavior and has been widely recognized in the LST measurement. The effect of thermal radiation directionality (TRD) shows that different view angles lead to different LST measurements [106]. Therefore, different thermal TDR models, including radiative transfer models, geometric models, hybrid models, 3-D models, etc., have been used to improve the LST measurement [107,108]. To better measure the emissivity and 3-D thermal anisotropy effects in the MODIS data, a series of building structure components, such as roofs, sunlit or shaded walls, and so on, should be integrated into the TRD models (e.g., CoMSTIR, SUM) to simulate the thermal infrared radiation of 3-D urban buildings [109,110]. MODIS observations are also influenced by the anisotropy effect [111,112], which is most evident for images with view angles of around 55° over complex terrains. The data pre-processing procedure in this study was not corrected for this effect; therefore, the summer mean LST was calculated for each mega-region to mitigate its impact in order to ensure that the anisotropic distribution of thermal radiance was on average well represented.

Li et al. [55] demonstrated that the LST-UDI correlation depends on the zoning scale of built-up areas or the spatial resolution of satellite images. Therefore, the scale issue should be considered when analyzing correlations between LST and landscape features [113]. In addition, the UDI definition in this study was based on the density of ISA instead of the urban-rural boundary, and thus cannot further demonstrate the interaction between the urbanization process and the urban heat hotspot. The impacts of other influencing factors of urban heat phenomena (e.g., increased emissions of anthropogenic heat, changes to urban geometry and landscape configuration, persistent high-pressure systems, moisture fluxes, and seasonal climate variability) are not within the scope of this study, but should be discussed in further studies.

Last but not the least, though there is a strong relationship between LST and NSAT, the two temperatures have different physical meanings, magnitudes, measurement techniques, diurnal phases, and responses to atmospheric conditions [112,114]. In this study we directly applied indicators developed for NSAT to measure the LST-based UHW phenomena. The differences of UHW characteristics, especially the coupling effect of LST and NSAT under global climate change, should be taken into consideration in the future.

6. Conclusions

Rapid urbanization has significantly reshaped the urban thermal environment, but conventional studies using in situ air temperature observations are often constrained by the limited coverage of the observation network. In this study, the MODIS land surface temperature (LST) product was used as a proxy for air temperature to depict the urban heat island (UHI) and urban heat wave (UHW) in nine mega-regions globally between 2003 and 2020. A daily window was used to calculate the combined hot days and tropical nights (CHT) based on daily minimum and maximum temperature to depict the frequency of absolute threshold-based UHW events, and a 3-day window was adopted to depict the percentile-based heat wave number (HWN) in the built-up and Non-built-up areas, respectively. On that basis, the correlation between heat wave magnitude or intensity and built-up density was further analyzed using impervious surface area data.

From 2003 to 2020, an increase in the mean LST as well as LST spatial variability was observed in all the mega-regions. Overall, a 1 °C increase in UHI intensity corresponds to a 2.04% to 92.15% increase in CHT, as well as a 3.30% to 33.67% increase in HWN. This suggests that under the combined influence of increasing LST as well as LST variability, the high temperature extremes represented by the increase in UHI magnitude and UHW frequency would occur more frequently under global warming.

The urban built-up density has a strong impact on urban heat phenomena. Overall, a 10% increase in the urban built-up density corresponded to a 0.20 °C to 0.95 °C increase in LST, a 0.59% to 7.17% increase in CHT, as well as a 0.08% to 0.95% increase in HWN. Over half of the impervious surface areas in the mega-regions have experienced at least one heat wave event. The increasing rates of CHT were higher than the expansion rate of built-up areas in six mega-regions, particularly in Boston, Tokyo, São Paulo, and Mexico City.

Percolation theory was adopted to delimit the hotspots of urban thermal environment. Using LST as the indicator, the hotspot magnitude ranged from 97.32% to 88.13%, while the UHW-based hotspot magnitude ranged from 86.34% to 1.18%. The increasing magnitude of UHI/UHW hotspots reveals that the urban expansion had—and will continue to have—detectable influence on the thermal environment at the mega-region scale. Therefore, the spatial correlation between the urban impervious surface and the urban thermal environment should be considered in future research to explore the mechanistic drivers of heat waves.

In this study, LST was utilized as a proxy for air temperature. Though they are closely correlated and consistent under most conditions, their differences should be considered in the interpretation of the results. The urban heat environment is also affected by factors such as land–atmosphere coupling, soil wetness condition, and the enhanced greenhouse effect, which were not analyzed in this study.