The Impact of Climate Change on Urban Thermal Environment Dynamics

: The human population is increasing. The ongoing urbanization process, in conjunction with climate change, is causing larger environmental footprints. Consequently, quality of life in urban systems worldwide is under immense pressure. Here, the seasonal characteristics of Maribor’s urban thermal environment were studied from the perspectives of surface urban heat island (SUHI) and urban heat island (UHI) A remote sensing thermal imagery time series and in-situ measurements (stationary and mobile) were combined with select geospatial predictor variables to model this atmospheric phenomenon in its most intensive season (summer). Finally, CMIP6 climate change scenarios and models were considered, to predict future UHI intensity. Results indicate that Maribor’s UHI intensity maximum shifted from winter to spring and summer. The implemented generalized additive model (GAM) underestimates UHI intensity in some built-up parts of the study area and overestimates UHI intensity in green vegetated areas. However, by the end of the century, UHI magnitude could increase by more than 60% in the southern industrial part of the city. Such studies are of particular concern, in regards to the increasing frequency of heat waves due to climate change, which further increases the (already present) heat stress in cities across the globe.


Introduction
Nowadays, cities are the dominant habitats of the human species. In 2000, about 47% of people lived in urban systems [1]. This proportion increased to 55% by 2018 and could reach 60% by 2030 [1]. By then, it is expected that 44% of the world's population will live in cities of at least half a million inhabitants [1]. With urban area expansion, artificial materials (concrete, asphalt, tiles, metals etc.) replace natural vegetation [2]. Urban construction thus causes a change in surface albedo, emissivity, and heat capacity [3]. As a result, a net influx of heat into the urban atmosphere is detectable [4], which is further potentiated by anthropogenic activities [5] in which heat is released (heating, industry, and transport) [6]. This evident positive air temperature anomaly between urban and rural areas is the urban heat island (UHI) [7]. UHI is considered one of the most problematic consequences of rapid urbanization [7]. It can cause health issues [5] and decrease productivity [8], especially during increasingly frequent, climate change-related heat wave events [9].
To measure and evaluate UHI intensity, two methodological approaches prevail. Larger metropolitan areas and mid-size cities require thermal remote sensing data and, thus, cloudless meteorological conditions, in order to study the surface urban heat island Maribor lies at the conjunction of the pre-Alpine and Pannonian areas of northeastern Slovenia (Figure 1a), along the river Drava. The city, with 95.767 inhabitants today [21], flourished in the middle of the 19th century, when the first industrial plants were built in connection with the construction of the railway network program between Vienna (AUT) and Trieste (ITA). Soon, other industries began to develop, but experienced a decline after the breakup of Yugoslavia in 1991. Today, Maribor is attempting to create its identity in the service economic sector and as a university center for northeastern Slovenia. The city has experienced accelerated population growth and spatial development since the 1960s, mainly owing to intensive rural immigration. The built-up part spread to the south, west, and east, occupying former farmland. Today, several parks (Mestni park, Trg Borisa Kidriča, Slomškov trg, Magdalenski park) and two forest areas (Stražun and Betnavski gozd) are still preserved (Figure 1b). Moreover, the river Drava divides Maribor into two parts, with a riverbed up to 160 m wide; it is thus an important local climate modifier [22].

Seasonal SUHI Characteristics
Surface thermal properties are among the main contributors to the UHI phenomenon. In order to capture seasonal Land Surface Temperature (LST) variability in the study area, 22 LANDSAT 8 cloud free satellite images (5 in winter (December, January, February), spring (March, April, May), and summer (June, July, August), and 7 in autumn (September, October, November)) between 2013 and 2020 were obtained from the EarthExplorer web platform [23] (path = 190, row = 28; scene center time between 10 and 10:25 am GTM). Both thermal channels (band 10 and 11) were initially converted to at-satellite brightness temperature (in • C) by applying TerrSet2020's Landsat algorithm [24] and were later averaged. Seasonal mean values of the resulting LST variable were calculated in the next Figure 1. The geographic position of the study area (a); stationary measurements (blue dots, 1 = the medieval city center (Industrial, commercial, public, military, and private units); 2 = the Melje industrial zone (Industrial, commercial, public, military, and private units); 3 = green urban area Mestni park (Green urban areas); 4 = the Tezno industrial zone (Industrial, commercial, public, military, and private units); 5 = Stražun forest (Forests); 6 = rural background (Arable land)), and the mobile measurement route (red lines) (b).

Seasonal SUHI Characteristics
Surface thermal properties are among the main contributors to the UHI phenomenon. In order to capture seasonal Land Surface Temperature (LST) variability in the study area, 22 LANDSAT 8 cloud free satellite images (5 in winter (December, January, February), spring (March, April, May), and summer (June, July, August), and 7 in autumn (September, October, November)) between 2013 and 2020 were obtained from the EarthExplorer web platform [23] (path = 190, row = 28; scene center time between 10 and 10:25 am GTM). Both thermal channels (band 10 and 11) were initially converted to at-satellite brightness temperature (in °C) by applying TerrSet2020's Landsat algorithm [24] and were later averaged. Seasonal mean values of the resulting LST variable were calculated in the next step. Finally, we were able to estimate SUHI intensity according to the land use/land cover categories of the Urban Atlas (UA) 2018 database, provided by the Copernicus web platform [25]. The category Pastures (code 23000) served as the reference.

Seasonal UHI Characteristics
These two datasets (SUHI intensity and UA) were the basis for placing 6 HOBO MX2302A temperature and moisture loggers (accuracy = +/− 0.2 °C and +/− 2.5% RH), in Davis radiation shields, in different local climate zones [26] of the city. We covered 4 UA categories (arable land (1), forests (1), green urban areas (1) and industrial, commercial, public, military, and private units (3, considering the build-up area gradient)) ( Figure 1b). All sensors were installed at 3 m height on an azimuth of 180° after preliminary testing of possible differences concerning the standard meteorological measurement height (2 m). A key reason for choosing a non-standard height was protection against vandalism. Air temperature and relative humidity data were recorded at half-hour intervals. The measurements began on 10 September, 2019. Thus, this study considers 659 measurement days and all four seasons (spring, summer, autumn, and winter). Fixed loggers provided atpoint thermal boundary layer atmospheric conditions in different parts (local climate zones) of the city, but we were still lacking spatial information about the UHI phenomenon. To bridge this gap, mobile air temperature measurements were performed in summer 2020 and 2021, with an IoT MF-300 temperature datalogger (GPS position accuracy = +/− 2 m; GPS velocity accuracy = 0.1 m/s; GPS center frequency = 1575 +/− 3 MHz; temperature accuracy = +/− 0.3 °C; temperature resolution = 0.1 °C; acquired at www.pileus.si, accessed on 1 September, 2019). Its -165 dBm tracking sensitivity extends positioning cov-

Seasonal UHI Characteristics
These two datasets (SUHI intensity and UA) were the basis for placing 6 HOBO MX2302A temperature and moisture loggers (accuracy = +/− 0.2 • C and +/− 2.5% RH), in Davis radiation shields, in different local climate zones [26] of the city. We covered 4 UA categories (arable land (1), forests (1), green urban areas (1) and industrial, commercial, public, military, and private units (3, considering the build-up area gradient)) ( Figure 1b). All sensors were installed at 3 m height on an azimuth of 180 • after preliminary testing of possible differences concerning the standard meteorological measurement height (2 m). A key reason for choosing a non-standard height was protection against vandalism. Air temperature and relative humidity data were recorded at half-hour intervals. The measurements began on 10 September, 2019. Thus, this study considers 659 measurement days and all four seasons (spring, summer, autumn, and winter). Fixed loggers provided at-point thermal boundary layer atmospheric conditions in different parts (local climate zones) of the city, but we were still lacking spatial information about the UHI phenomenon. To bridge this gap, mobile air temperature measurements were performed in summer 2020 and 2021, with an IoT MF-300 temperature datalogger (GPS position accuracy = +/− 2 m; GPS velocity accuracy = 0.1 m/s; GPS center frequency = 1575 +/− 3 MHz; temperature accuracy = +/− 0.3 • C; temperature resolution = 0.1 • C; acquired at www.pileus.si, accessed on 1 September, 2019). Its -165 dBm tracking sensitivity extends positioning coverage into places such as urban canyons and dense foliage environments, where GPS usage was not previously possible. The datalogger's air temperature sensor was installed on a vehicle radio antenna (50 cm high) and was additionally protected with a small Davis radiation shield. Data were recorded each second (Figure 1b). Because of the road network structure, at certain points, such as road traffic lights, air temperature measurements were repeated. To filter out this noise (repeated air temperature records), a Python-based software was developed, which considered the exact time of the measurement and its direct spatial vicinity. Thus, other data errors were removed as well. Moreover, we measured air temperature differences, based on our stationary measurements approximately at the time when they reached the maximum value above urban structure (16:00 local time; started at 15:30, and finished at 17:00). Consequently, the mobile measurement time interval did not disrupt the results since differences in temperature records between the start and end of the measurement were statistically insignificant and, thus, ignored in further analysis. Because only temperature differences were analyzed, the ventilation effect (and traffic load) due to car movement (average speed = 40 km/h) was also disregarded. Altogether, 22 summer measurement days (90.444 air temperature records) were analyzed in different weather conditions/types (anti-cyclonic, cyclonic, and advective).

UHI Modeling
These mobile air temperature difference measurements enabled the development of the dependent variable summer UHI intensity (UHIi) in Maribor. To do so, a 300 m grid was designed in the ArcGIS environment [27] and laid over the entire study area. Next, the average air temperature for each quadrate/polygon covering the measuring route was calculated. Quadrates/polygons with less than 30 air temperature records were ignored. The Zonal statistics algorithm then enabled air temperature data extrapolation to UA categories. UHIi was calculated by subtracting the mean temperature value of the Pastures category (code 23000) from all other categories (as in the case of SUHI intensity). Finally, 10 random points per UA category (190 in total) were generated to operate with the point data, properly fit predictive models, and evaluate potential climate change impacts on Maribor UHI intensity.
To explain and model the detected UHIi geospatial pattern, the predictor variables per UA category were then produced: summer SUHI intensity (SUHIi), southness (= sinAspect) (sness), summer NDVI (NDVI), distance to water bodies (d2wb), building volume (bv), high vegetation volume (hvegv), and summer daily mean air temperature (June, July, August) (tas_summer). The predictor SUHIi was designed in the initial step as the product of LADSAT-based LST in reference to the Pastures category of the UA. The explanatory variable sness (a linearized form of the aspect variable where azimuth values are transformed into an index reaching values between −1 (northern slopes) and 1 (southern slopes)) was calculated in the ArcGIS environment from a LIDAR-based digital terrain model with a 1 m horizontal resolution, provided by the Slovenian Environmental Agency, which operates under the Ministry of the Environment and Spatial Planning [28]. Mean NDVI values per UA category were obtained by processing (Band 8-Band 4/Band8 + Band4) a Sentinel 2A satellite image (8 August, 2020) gathered from the Copernicus Open Access Hub [29]. Water bodies were extracted from the national land use database, which is freely available on the web platform owned by the Ministry of Agriculture, Forestry, and Food [30]. To transform this categorical variable, the Euclidean distance algorithm was applied in the ArcGIS environment [27]. Two databases were used to prepare the bv predictor. The above-mentioned LIDAR point cloud and the National register of buildings (a digital vector database), provided by the Ministry of Infrastructure and Spatial Planning [31]. A normalized digital surface model (nDSM) was calculated by considering the difference between first and ground returns. The National register of buildings was then used to calculate building volume per UA category. A similar approach was used to calculate high vegetation (tree) volume (hvegv). Vegetation cover was determined with NDVI and additionally validated with LIDAR cloud data classification (low, medium, or high vegetation). The tas summer predictor was developed for four time windows (current (1981-2011), 2011-2040, 2041-2070 and 2071-2100), two climate change scenarios (ssp370 (level of forcing = 7.0 W/m 2 in 2100) and ssp585 (level of forcing = 8.5 W/m 2 in 2100)) and five global climate models (GCMs) (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, UKESM1-0-LL), as provided by the CHELSA V2-CMIP6 database (spatial resolution = 30 arc seconds) [32]. To simplify modeling procedures, we used the tas_summer ensemble mean value from five models per future time window.
These predictors were finally checked for possible multicollinearity (Pearson's correlation coefficient > 0.5 or < −0.5) and later used to model Maribor's summer UHI intensity in current and future climate conditions. A generalized additive model (GAM) [33][34][35][36] was calibrated in the R statistical environment [37] by applying the mgcv package [38].
GAMs provide a structure for generalizing a general linear model by allowing additivity of non-linear functions of the variables [39]. The advantage of GAM is to limit the error in prediction of a dependent variable from various distributions by assessing unspecific functions, which are connected by means of a link function with the dependent variable [40]. Thus, our GAM model adequacy was tested with a generalized linear model (GLM) algorithm with the same predictors and designed with the Rcmdr package [41]. However, to additionally assure GAM model quality, its standardized residuals (overand under-predictions) were tested for normal distribution and spatial autocorrelation (Moran's I) in the ArcGIS environment [27]. After properly specifying the GAM model, the predict function within the mgcv package enabled future UHIi estimation by considering GCMs-based tas_summer values in time windows 2011-2040, 2041-2070, and 2071-2100, and emission scenarios ssp370 and ssp585.

Modeling Limitations
Since we sought to evaluate how future thermal conditions could influence the existing summer UHI pattern in Maribor, only the tas_summer variable was changed in the future climate change scenarios being considered. Of course, other UHI predictors will substantially change in the future and could alter the existing spatial pattern either positively or negatively. Higher air temperatures could potentially even decrease UHI intensity, especially in colder conditions (winter, spring, autumn). However, land use [42][43][44] and surface spectral change studies [45], as well as demographic trends [44,46] provided evidence that Maribor's urban structure will not be particularly dynamic in the future. This fact excuses our consideration of all remaining UHI intensity predictors (bv, d2wb, NDVI, sness, SUHIi, hvegv) as temporarily stable. Moreover, we decided to use such an approach to illustrate possible outcomes in the current city structure if no climate change adaptation and mitigation actions are implemented. In order to improve the climate change-related UHI pattern change estimation, city development plans should be integrated. This information could be used to prepare future estimates for several of these UHI predictors (i.e., bv, d2wb, NDVI, and hvegv) and, thus, produce even more accurate results.

The Seasonal SUHI Intensity Footprint
The spatiotemporal analysis of the LST variable revealed that Maribor's SUHI intensity pattern was most homogeneous and least pronounced in winter ( Figure 2). The Tezno industrial zone in the SE part of the city and the southern-exposed surroundings occupied by vineyards (Mariborske gorice) on the northern built-up edge, stand out from the reference areas. Temperature differences between other land use categories were smaller ( Figure 3). In spring, the SUHI intensity footprint increases and diversifies. Positive surface temperature anomalies were evident in the industrial zones of Tezno (SE part) and Melje (E part), in Studenci (W part), and in the shopping center area with associated parking lots. Negative deviations were observed along the river Drava and the corresponding hydroelectric canal, and in all closed forest areas (Stražun (E part), Betnavski gozd (S part), and Limbuška Dobrava (W part)). All of the smaller hot spots that had been identified in the spring were amalgamated in the summer. SUHI intensity increased up to 12 • C. All three industrial areas and the medieval center generated high levels of summer heat stress. Fortunately, the river Drava and other green urban areas create a large surface temperature contrast and, thus, break down the otherwise uniform SUHI footprint. In autumn, SUHI intensity decreases and simplifies. The spatial pattern is similar to that formed in spring.
and Limbuška Dobrava (W part)). All of the smaller hot spots that had been identified in the spring were amalgamated in the summer. SUHI intensity increased up to 12 °C. All three industrial areas and the medieval center generated high levels of summer heat stress. Fortunately, the river Drava and other green urban areas create a large surface temperature contrast and, thus, break down the otherwise uniform SUHI footprint. In autumn, SUHI intensity decreases and simplifies. The spatial pattern is similar to that formed in spring.

SUHI Intensity and Land Use
In the summer, the highest surface temperature differences compared to the reference UA category (pastures) occurred on continuous urban fabric areas (sealing level (SL) > 80%) (by 4.0 °C), on industrial, commercial, public, and military surfaces (by 3.7 °C), on railways, and associated land (by 3.3 °C) (Figure 3). Surface temperatures were lower on water surfaces (by 2.6 °C), in forest areas (by 1.8 °C), and for isolated structures (by 0.2 °C). A similar, but less intensive, pattern was observed in the spring and autumn. In winter, the highest positive deviations in LST occurred in the category permanent crops (represented mainly by vineyards) (by 0.5 °C). These are located on the steep southern slopes

SUHI Intensity and Land Use
In the summer, the highest surface temperature differences compared to the reference UA category (pastures) occurred on continuous urban fabric areas (sealing level (SL) > 80%) (by 4.0 • C), on industrial, commercial, public, and military surfaces (by 3.7 • C), on railways, and associated land (by 3.3 • C) (Figure 3). Surface temperatures were lower on water surfaces (by 2.6 • C), in forest areas (by 1.8 • C), and for isolated structures (by 0.2 • C). A similar, but less intensive, pattern was observed in the spring and autumn. In winter, the highest positive deviations in LST occurred in the category permanent crops (represented mainly by vineyards) (by 0.5 • C). These are located on the steep southern slopes in the hilly surroundings of the city (Mariborske gorice). In contrast, forest areas, water bodies, arable land, and green urban areas were significantly cooler than pastures in winter.

UHI Seasonal Dynamics
Remotely sensed data revealed the seasonal surface thermal pattern of the study area. Based on these findings, atmospheric thermal conditions were captured in different local climate zones (four UA categories (arable land, forests, green urban areas and industrial, commercial, public, military, and private units)). Figure 4 summarizes the seasonal development of the UHI phenomenon in Maribor. The two-way ANOVA analysis and the corresponding Tukey post-hoc test (p < α; α = 0.05) additionally confirmed significant differences in UHI intensity between seasons and UA categories (local climate zones). However, the largest differences in daily mean air temperature were evident in summer. The Stražun forest was clearly cooler than the reference surface (arable land in the city background), followed by the Mestni park green urban area, which was occasionally even hotter. Instead, all three built-up surfaces (the Melje industrial zone, the medieval city center, and the Tezno industrial zone) were significantly hotter than the rural surroundings. The same pattern was observed throughout the year, in spring and summer, with UHI intensity extremes occurring mainly in spring. the Tezno industrial zone) were significantly hotter than the rural surroundings. The same pattern was observed throughout the year, in spring and summer, with UHI intensity extremes occurring mainly in spring.

The Summer UHI Model
Both the perspectives on Maribor's thermal environment (SUHI and UHI) indicated that citizens were exposed to extreme heat stress, particularly in summer. Thus, the focus of the research was turned to the summer season UHI. A GAM predictive model was designed by considering several explanatory variables known for their impacts (linear or non-linear) on the UHI phenomenon. Table 1 shows that all the given predictors met the basic criteria, since correlation coefficients were within the −0.5 and +0.5 margin. How-

The Summer UHI Model
Both the perspectives on Maribor's thermal environment (SUHI and UHI) indicated that citizens were exposed to extreme heat stress, particularly in summer. Thus, the focus of the research was turned to the summer season UHI. A GAM predictive model was designed by considering several explanatory variables known for their impacts (linear or non-linear) on the UHI phenomenon. Table 1 shows that all the given predictors met the basic criteria, since correlation coefficients were within the −0.5 and +0.5 margin. However, low variance inflation factors (VIF < 3) additionally excused further use of these predictors in the modeling procedure. The GAM approach fitted with seven predictor variables yielded 82.4% of explained deviance ( Table 2). All predictors, except tas_summer, had a statistically significant impact (p < α; α = 0.05) on the summer UHIi pattern. However, we kept the tas_summer predictor in the model, since it described current (and in predict mode, future) summer daily mean air temperatures in the study area. Four predictors (SUHIi, sness, NDVI, and d2wb) were non-linearly related to the dependent variable (edf > 1). The scale estimate (residual standard error squared = 0.013) indicated good predictive power for the designed UHIi GAM model. Nonetheless, we additionally tested its quality by comparing it against a GLM approach, which assumes that all explanatory variables are linearly related with the given dependent variable ( Table 3). The analysis of deviance, based on χ 2 statistics, revealed that Maribor's UHIi can be better modeled by applying a GAM (p < α; α = 0.05). However, before reaching any conclusions, standardized residuals were tested for possible spatial autocorrelation. The insignificant Moran's I index (p > α; α = 0.05) proved that model over-and under-predictions were normally distributed and free of spatial autocorrelation ( Figure 5). In other words, the model was properly specified without missing predictor variables. Moreover, Figure 5 indicates that major over-predictions are linked to the rural surroundings of the city, whereas under-predictions relate to the urban structure.  Finally, the developed GAM, the predict function in the mgcv package and the prepared climate scenario tas_summer predictor enabled future UHIi estimation in Maribor.  Finally, the developed GAM, the predict function in the mgcv package and the prepared climate scenario tas_summer predictor enabled future UHIi estimation in Maribor.   (current (1981-201 2011-2040, 2041-2070, and 2071-2100) indicate a gradual increase in UHIi. Some possi improvements in summer UHIi are predicted only in the first half of the century (main under ssp370 conditions). The second half could potentially bring up to a 60% (> 0.5 ° UHIi increase, in particular in the southern industrial part of the city. However, diff ences between emission scenarios are evident ( Figure 6).

Discussion
Ensuring high quality of life in urban systems is becoming increasingly difficult. In the last half of the century, the share of urban population has risen from 34% to 54%, and it is expected to reach 66% by 2050 [47]. In addition, cities worldwide are, owing to climate change, already exposed to higher mean air temperatures, changed precipitation patterns, and more frequent extreme weather events [48,49]). In NE Slovenia, where Maribor is located, air temperatures rose between 0.33 and 0.44 • C per decade (period 1961-2016) and in the summer months, between 0.44 and 0.59 • C per decade [50,51]). There, heat wave frequency has more than tripled since the 1960s [22]. Moreover, heat wave events are lasting longer (8.8 days more since 1990) and reaching higher maximum temperatures (up to 40.6 in 2013). However, urban warming is a manifestation of the direct (sensible heat) and indirect (land-cover transformation, climate change) alteration of the energy budget in the urban boundary layer [4,52,53]. These facts lead to the conclusion that a consequent intensification in the UHI effect is inevitable, in large agglomerations or small cities/towns, if no preventive or adaptation actions are implemented. From that perspective, studies like this, which provide applicative spatial prediction models, and which could roughly estimate the potential change in UHI intensity owing to direct or indirect urban system alterations, are of high importance. However, the annual or seasonal dynamics of the UHI phenomenon in small cities/towns under climate change scenarios is, in general, poorly studied. We sought to bridge this gap and, thus, broke down Maribor's thermal environment into individual entities that cause the current spatial UHI pattern, and will potentially co-create its thermal footprint in the second half of the 21st century. A time series of multispectral satellite imagery, along with stationary and mobile air temperature measurements, provided the first important discovery: the UHI intensity maximum in Maribor shifted, owing to climate change, from winter to spring and summer. Several urban climate scientists across the globe [4]) agreed that the UHI phenomenon in mid-latitude cities/towns is usually more pronounced in winter, moments before sunrise (longer nights provide more time for cooling and intensive house heating causes higher maximum air temperatures). The applied GAM summer UHI intensity model based on seven spatial predictor variables yielded good results. It should be emphasized that we tried to model Maribor's summer UHI intensity with a multiscale geographically weighted regression (MGWR) [54], alongside the above-mentioned GLM. It yielded slightly better results (adjusted R2 = 0.84 and smaller over-and under-predictions), but we decided to calibrate a GAM because local regression models are still case-specific and cannot easily be transferred elsewhere in the geographical space [15,55]. Thus, our second important finding is that major over-predictions are linked to the rural surroundings of the city, whereas under-predictions exist in relation to the urban structure. In other words, the GAM model underestimates UHI intensity in some built-up parts of the study area and overestimates UHI intensity in green vegetated areas. Future UHI intensity scenarios could thus potentially be even worse. However, it should be emphasized that our predictions are made under the assumption that the city structure (predictor variables bv, d2wb, NDVI, sness, SUHIi, hvegv) remains temporarily stable. In other words, the uncertainty level of our predictions increases as soon as land use changes are applied. Other important summer UHI intensity modifiers that are not under consideration could decrease the proposed GAM model predictive power as well. For instance, a changed summer rainfall pattern could affect water bodies, vegetation (NDVI), and wind speed. As already mentioned, rising air temperatures alone could also decrease UHI intensity, especially in colder conditions (winter, spring, autumn). Studies about future bioclimatic conditions in the NE part of Slovenia (where Maribor is located) indicated that this region could potentially face hotter and dryer future climate conditions [56,57]. From that perspective, our results remain highly applicable despite the outlined weaknesses. Nonetheless, the third finding is that our methodology can be transferred to other cities across the globe, since it exploits statistical relations over all the data in the region of interest, and can therefore be applied even beyond the geographical limits of the measurements in cases where the physical settings remain similar. Finally, the new CMIP6 high-resolution climate change scenarios [32] indicate that Maribor could soon face a reasonable increase in summer UHI intensity, in particular in the southern industrial zone. By the end of the century, summer UHI magnitude in this area could increase by more than 60%. From that perspective, the last important message of this study is as follows: climate change preventive actions in urban systems are urgently needed to mitigate high levels of heat stress and assure a quality living environment for all social groups of citizens. Predictions for Maribor, which has several green urban areas and a large water body (the river Drava), are not optimistic. This is even more alarming for other Slovenian urban settlements, which are for historical reasons [58], much more concentrated.
Climate change and the UHI phenomenon are linked through a two-way interaction. First, global warming will increase UHI intensity in urban systems. Second, cooling strategies to reduce UHI intensity can help communities adapt to the impact of climate change, while also lowering the greenhouse gas emissions that cause climate change [59]. Most efforts to cool UHIs produce many benefits, including lower temperatures, greenhouse gases, air pollution, harmful health impacts, and electricity demand [60][61][62]. Efforts to reduce elevated temperatures in urban systems thus also help to address climate change and improve air quality. Moreover, such measures can help citizens become more resilient to some of the damaging impacts of modern climate change. Properly located green areas and water bodies represent good examples for urban cooling [63][64][65][66]. In this way, the need for fossil fuel-based transport (which is still dominant) to cooler parts of the city is reduced. In order to decrease the record-breaking CO 2 emissions (417 ppm, accessed on 17 June, 2021) [67] and somehow mitigate the impact of climate change at the global level, several multifaceted mitigation and adaptation measures must be taken as soon as possible. Starting with UHI mitigation strategies in densely populated areas is certainly a prime target, since countermeasures would, thus, be implemented precisely at the source of the greenhouse gas emissions.