Climate Justice in the City: Mapping Heat-Related Risk for Climate Change Mitigation of the Urban and Peri-Urban Area of Padua (Italy)

: The mitigation of urban heat islands (UHIs) is crucial for promoting the sustainable development of urban areas. Geographic information systems (GISs) together with satellite-derived data are powerful tools for investigating the spatiotemporal distribution of UHIs. Depending on the availability of data and the geographic scale of the analysis, different methodologies can be adopted. Here, we show a complete open source GIS-based methodology based on satellite-driven data for investigating and mapping the impact of the UHI on the heat-related elderly risk ( HERI ) in the Functional Urban Area of Padua. Thermal anomalies in the territory were mapped by modelling satellite data from Sentinel-3. After a socio-demographic analysis, the HERI was mapped according to ﬁve levels of risk. The highest vulnerability levels were localised within the urban area and in three municipalities near Padua, which represent about 20% of the entire territory investigated. In these municipalities, a percentage of elderly people over 20%, a thermal anomaly over 2.4 ◦ C, and a HERI over 0.65 were found. Based on these outputs, it is possible to deﬁne nature-based solutions for reducing the UHI phenomenon and promote a sustainable development of cities. Stakeholders can use the results of these investigations to deﬁne climate and environmental policies.


Introduction
The Intergovernmental Panel on Climate Change Sixth Assessment Report [1] states that in the next 20 years, the global temperature will increase by 1.5 • C or more.Recordbreaking temperatures will be more frequent, intense and last longer due to climate change [2].As climate change proceeds, urban heat-related health risks are likely to intensify over the next decade.In fact, urban heat islands (UHIs) have a significant impact on health by increasing mortality [3].Recently, Europe was identified as a heat wave (HW) hotspot, and it is speculated that temperatures may rise three to four times faster than in the rest of the northern mid-latitudes [4].
The UHI phenomenon is strongly influenced by the land cover/land use and the builtup environment.The reduction in vegetated areas and urban expansion promote an increase in the land surface temperature (LST), thus exacerbating the UHI phenomenon [5,6].
Since the number of people in urban areas tends to increase, it is essential to promote climate-resilient development to improve liveability and reduce risks related to atmospheric conditions, such as climate change and UHI.In applying adaptation policies and strategies to mitigate the impact of extreme events, the most degraded areas are often penalised.
To date, few studies have documented obstacles to climate justice or assessed how emerging adaptation plans impact marginalised groups.To address the issue of climate justice, participation in adaptation planning should be broadened, with a focus on fastgrowing/low-budget cities, by integrating justice into urban design processes [7].
Defining methodologies for reducing the LST and spatially assessing UHI hotspots could help urban planners and decision makers develop climate-resilient policies grounded in nature-based solutions (NBSs) as well as participative and inclusive processes.

Research Background
Various methods and techniques for assessing UHI can be found in the literature.The main criteria for defining the methodology largely depend on the geographical scale of the analysis and the availability and quality of data.
The following subsections describe (i) the methods used to assess and quantify the UHI phenomenon, (ii) the strategies for UHI mitigation and assessing their impacts on the territory (iii) and a description of some representative study cases.
In general, it is necessary to develop an optimised analysis method for UHI estimation [8].Currently, there are a lack of protocols for establishing UHI classification criteria and for determining data collection methods that can efficiently estimate and analyse UHIs [9].

Urban Heat Island (UHI) Assessment
The most commonly used methods for assessing UHIs are integrated geographic information system (GIS) and remote-sensing techniques based on processing satellite images [10,11].
Regarding the input data (i.e., satellite images), Landsat satellite images are most often used to estimate the LST and land cover; thermal bands are used (i.e., TIRS1 and TIRS2 of Landsat 8) at a spatial resolution of 100 m.However, their application may be limited due to satellite revisit time and, therefore, cloud cover.Sentinel satellite images can address this issue, as they provide high-temporal-frequency data, but they have a spatial resolution of 1 km [12].More specifically, the Sea and Land Surface Temperature Radiometer (SLSTR) onboard Sentinel-3 satellites provides high-quality observations that can be used to estimate the LST.Sentinel-3 images have been used in various studies, despite their low spatial resolution.For example, Yang et al. [13] defined a series of algorithms to estimate the LST at six sites in northwest China, using Sentinel-3 images and land cover data with a spatial resolution of 1 km 2 .In another work, Landsat 8 and Sentinel-3 images were used in combination to calculate the LST through algorithms.The investigation was carried out in the metropolitan area of Granada, which has a territorial surface area of 880,000 km 2 .In both cases, reliable results were obtained [14].Other works [15,16] have also confirmed that Sentinel-3 provides reliable results.In general, Sentinel-3 LST values are slightly higher than LST values measured in situ.Therefore, based on the geographical and temporal scale on which analyses are performed, certain data might be more suitable than others.
Different GIS software programs are currently in use.ESRI ArcGIS ® is one of the most popular for spatial data processing and modelling despite being a paid, proprietary program.Open source software programs, such as QGIS, have shown similar performance and are increasingly popular [17].
Currently, the main methods for assessing UHI phenomenon are historical data from ground weather stations, computer modelling and simulations and remote-sensing technologies.Correlations between UHI-related variables and UHI phenomenon using remotesensing data are the most commonly used technique [11].Moreover, by carrying out sensitivity analyses, it is possible to identify the parameters that influence the development of UHIs.These parameters include the albedo, the presence of vegetation, the population density and the urban morphology, which is described using street widths, canyon orientation and building heights [18,19].

Mitigation Strategies and Their Impacts
Green infrastructures (GIs) can be used to mitigate the effect of UHIs.Through evapotranspiration and shading, GIs improve outdoor thermal comfort.For example, Emmanuel and Loconsole [20] predicted that a 20% global increase in green cover could reduce surface temperatures by 2 • C by 2050.
More specifically, there are different GI strategies that are implemented in spatially based urban planning, including the use of trees (the most used GI strategy), grass, shrubs, green roofs and walls and parks [21].A tree cover of at least one-third of the total surface area reduces the temperature by about 1 • C [22], and the use of trees can reduce the temperature, on average, from a minimum of 0.5 to a maximum of 1.5 • C [21,23,24].The combined use of trees and grass can reduce the external air temperature up to 2-2.3 • C [25,26].Meanwhile, the use of shrubs generates a moderate, sometimes minimal reduction in temperature.The impact of green roofs on temperatures is significantly lower than that of trees, in some cases even nil, whereas green walls have a greater impact on reducing temperatures [22,27].Finally, green parks, have significantly lower temperatures than the urbanised surrounding area, with differences of 5 to 12 • C, and they can also reduce the temperatures of the surrounding buildings [28].
In general, the benefits of GIs on temperatures vary considerably depending on the geographical scale (city or building level), the extent of the greenery (shape and size of the park) and the type of vegetation [29].Further, to be effective and usable, mitigation strategies must be defined in conjunction with climate policies [30].However, they are usually defined separately, without interaction.

Case Studies
Almost all studies have used a remote-sensing approach, although some studies have used a more traditional approach.For example, Martinelli et al. [31] investigated the UHI in a Mediterranean coastal city, Bari, in Italy.As input data, the authors processed the urban canopy layer and the air temperature data recorded by four weather stations located in the metropolitan area of the city.
Conversely, the following studies used a remote-sensing approach and spatially based software to process data and investigate climate scenarios, such as GIS and ENVI-met © .Cocci Grifoni et al. analysed the impacts of climate change on the UHI in Ascoli Piceno City using remote sensing, spatial data, satellite images and computational fluid dynamics (CFD) simulations (i.e., ENVI-met © ) [32].Kikon et al. [33] used Landsat satellite images to assess the UHI in the city of Noida in India.In particular, a spatiotemporal analysis of LST was carried out, and correlations between LST and other land cover indicators, such as the Normalised Difference Vegetation Index (NDVI) and albedo, were found.In another work [34], Landsat 8 satellite images were used to study not only the LST and NDVI but also the Normalised Difference Built-up Index (NDBI).This method was applied to the Colombo metropolitan area in Sri Lanka.In a case study near Colombo in Kurunegala, NBSs for mitigating the UHI were defined using Landsat data series from 1996 to 2019 [35].A similar remote-sensing approach was applied to a Bogor case study to assess the LST distribution from 1990 to 2017 [36] and to a Guangzhou case study in China to investigate the spatial distribution of the LST and the anthropogenic heat flux from 2004 to 2020 [37].For these analyses, Landsat 5 and Landsat 8 satellite images were processed.Similarly, Landsat 8 images and land cover data were used in an analysis of Paris and Geneva [38].

Study Area: The Functional Urban Area (FUA) of Padua
The area under study is the urban and peri-urban area of Padua (NE of Italy), which include 31 municipalities, with a total population of 530,322 (data updated to 2020) and a territorial area of over 600 km 2 .We adopted a geographical scale, including whole functional urban area (FUA) of Padua into the spatial analyses.The European Commission identifies an FUA as an extension of a city's population density and work-related flows.In fact, an FUA consists of a densely populated city, in this case, Padua (209,730 inhabitants in 2020), and a surrounding area (commuting zone, 320,592 inhabitants in 2020).Figure 1 shows the municipalities included in the FUA and the reference official ground weather station (yellow dot in the municipality of Legnaro) used in this work to identify time windows of HWs (see Section 2.1).The local climate is humid subtropical, with cold winters and hot summers.This urban territory is one of the most affected in Italy by soil sealing due to urban expansion and infrastructure.In fact, according to an annual report by the National Institute for Environmental Protection [39], the city of Padua was ranked as the fifth largest city, with more than 100,000 inhabitants, with sealed surfaces.Recent studies that adopted a GIS-based approach estimated that about 50% of the municipal territory is currently sealed by urban infrastructure, jeopardising urban vegetation [40][41][42].
As has been widely documented [43][44][45][46][47], HWs, which affect urban areas that are extremely sealed by impervious surfaces, may exacerbate the impacts of UHIs on the area; therefore, it is crucial to map the phenomenon and identify mitigation policies.
In accordance with other studies [48][49][50], in Padua and the surrounding municipalities, there are uncomfortable conditions during the summer temperature peaks.Based on these studies, the air temperatures in these urbanised areas have increased by about 6 • C.

Research Objectives
This work proposes a replicable methodology to explore, assess and map UHIs and their impacts during heat wave periods.The general aim of this work is to investigate and assess the effects of HWs at the FUA geographical scale, using integrated GIS and remote-sensing techniques based on processing Sentinel-3 satellite images.Compared to other methodologies adopted, we combined the use of historical data based on ground weather stations together with remotely sensed data.This combination allowed us to firstly screen the timeframe for extreme heat events described by the heat wave magnitude index daily (HWMId), then to identify suitable scenes from Sentinel-3.Moreover, the use of GIS combined with open access satellite data such as Sentinel-3 allows to map and identify critical areas as a function of several variables, such as land cover/land use, surface temperatures, and socio-demographic data.In fact, data from remote sensing together with GIS modelling allow us to evaluate the temperature values, in this case, the LST, for the portion of the territory to be analysed (the resolution of the data varies according to the analysis scale, the availability of the data and the type of satellite).The specific objectives of this work are summarised as follows: UHI assessment by calculating the HW period and the heat wave magnitude index daily (HWMId) using a meteorological analysis from 2018 to 2021.II.Mapping the LST and NDVI and the percentage of green surfaces.III.Mapping urban thermal anomalies by calculating the temperature difference between the LST and the average LST in peri-urban green areas.IV.Mapping UHI hazard and vulnerability by modelling the heat-related elderly risk index (HERI).

Materials and Methods
The Materials and Methods Section is divided into five subsections.The first section describes the input data (i.e., weather data, satellite data, land use data, socio-demographic data) and the main tools used to process them.The other sections show the GIS-based methodology for evaluating the thermal anomalies and the impacts of the UHI on the HERI in the urban and peri-urban area of Padua.
Figure 2 shows the workflow of the GIS-based methodology; the main phases are: 1.
The collection of input data; 2.
The pre-processing of data using open source software, SNAP and QGIS; 3.
The processing phase, in which all the parameters necessary to assess the urban thermal anomalies and the HERI are calculated; 4.
In the last phase, the results are mapped in QGIS.
The proposed methodology is based on a complete open source and open data workflow, making it possible to be replicated and scaled up in other contexts.However, some constraints might exist related to the satellite scenes' relationship with the identified timeframe of extreme heat waves events.

Input and Processing Data
The input data and the tools used for geoprocessing and spatial modelling are free and open source.The weather data analyses were performed in Excel, the satellite data (Sentinel-3 images) were pre-processed in SNAP, and the geo-processing and the mapping of the results were developed in QGIS.Table 1 describes the spatial input data used in this work.For each layer, the source, the coordinate reference system (CRS) and a brief description are indicated.All data were re-projected into the reference system WSG 84/UTM zone 32N (EPSG: 32,632).Ground data refer to the 'Legnaro' weather station of the Regional Agency for Environmental Protection (ARPA), which is about 8 km from Padua, within the FUA.The hourly air temperature data refer to the last 30 years, from 1992 to 2021.In this work, we have considered the period 2018-2021 for the analyses, years for which Sentinel-3 satellite images are available.

Satellite Data
The use of Sentinel-3 was considered more appropriate than Landsat 8 for our research due to the geographical scale of the analyses and the availability of satellite images.The Sentinel-3 images were selected considering the HW time windows and the five warmest days of each year (see Section 2.2) from 1992 to 2021.Sentinel-3 images are available since 2016, so the data download refers to the period 2016-2021.Sentinel-3 satellite images downloaded for 2016 and 2017 were not used in this analysis due to satellite acquisition issues that made raster data not processible.Depending on the availability and quality of images, the analysis was conducted for the following days: 1 August 2018, 27 June 2019, 31 July 2020 and 14 August 2021 (24 June 2016 and 2 August 2017 were excluded due to data inaccuracies).
Sentinel-3 images were pre-processed in SNAP.They were re-projected using EPSG: 32,632 (WSG 84/UTM zone 32N) as a reference system and subsequently exported in tif format to be processed in QGIS.Since the study area corresponds to the FUA, the centroid of the FUA was calculated, and a 30 km buffer was created with respect to the centroid.The satellite images were extracted using the 30 km buffer as the clip area.The extracted parameters were NDVI and LST.LST values were converted from Kelvin to Celsius.Average NDVI and LST values were then calculated for used in subsequent analyses.In QGIS, the GRASS tool 'r.series' was used to calculate the average data of four scenes.
The LST values were mapped using Sentinel-3 data from four satellite scenes (summer period in the years 2018-2021).In particular, median LST data were calculated for each pixel with a resolution of 1 km × 1 km.The LST/NDVI values were mapped both using the original resolution of the Sentinel-3 image (Figures A1 and A3 in Appendix A), and the data were downscaled to a spatial resolution of 100 pixels for finer geovisualisation.To downscale data, the GRASS tool 'r.resamp.interp' in QGIS was used; in particular, data were processed according to the bilinear interpolation method.This process was carried out to improve the visualisation of the data and make it more readable.

Land Cover Data
The corine land cover (CLC) at 2018 (downloaded from the Copernicus database) was used to identify the vegetated surfaces.The latest version of land cover data is available for 2018; this is the third edition of CLC.The previous ones refer to the years 2007 and 2012 (updates are made every 5-6 years).The CLC we adopted was derived from high-resolution multispectral orthophotos (20 cm pixel size, 4 bands), at 1:10,000 nominal scale with a minimum mapping unit of 0.25 ha.Therefore, it was considered spatially and temporally appropriate for analyses carried out in the years from 2018 to 2021.The CLC database is in vector format, and the territorial surfaces are classified using a detail at the 5th level (the most detailed level available).The percentage of green surfaces was calculated with a spatial resolution of 1 km in order to compare these percentages with the LST and NDVI values processed by satellite images, which have the same spatial resolution.The first step was to create a grid (format: vector) with the same spatial resolution as the satellite images (1 km).A spatial intersection was performed between the grid and the CLC layer.Subsequently, it was possible to calculate the green area percentage for 1 km 2 units (cells).The vector data were converted to raster data, and the percentage of green surface was mapped for each pixel with a size of 1 km 2 .

Socio-Demographic Data
Socio-demographic information was collected from the statistical database of the Veneto Region at the municipal level (polygonal vector data), updated to 2020.
The population density (expressed in inhabitants per km 2 ) and the percentage of individuals over 65 years of age were calculated from this database.These data were then used to assess the HERI.

Urban Heat Island (UHI) Analysis
The UHI analysis was carried out by screening and identifying time periods affected by an HW and considering the five warmest days of the year.These days were identified using the 2 m air temperature recorded from the official weather station in Legnaro, 8 km from Padua.Starting with a statistical weather analysis of a 30-year time series (from 1992 to 2021) [51], the 2016-2021 period was selected because it is compatible with the availability of Sentinel-3 images.The most critical days from 2016 to 2021 were selected for the investigation.

Heat Wave (HW) Period
An HW is a period equal to or greater than three (method 1) or five (method 2) consecutive days with maximum air temperatures above the daily threshold for the reference period [52].The reference period considered is from 1992 to 2021.
The threshold is the 90th percentile of daily maxima, centred on a 31-day window.Therefore, for a specific day d, the threshold is the 90th percentile of the dataset A d defined by Equation (1): where: -U is the union of the datasets; -T y,i is the daily maximum temperature of the day i in the year y.

Heat Wave Magnitude Index Daily (HWMId)
The HWMId is an indicator used to calculate the severity of an HW [53].The HWMId is defined as the sum of the magnitude of the consecutive days comprising a HW (Equation ( 2)), with the daily magnitude calculated by Equation (3): with: where: -T d is the maximum daily temperature on day d of the HW; -T 30y25p is the 25th percentile of the annual maximum temperatures recorded from 1992 to 2021; -T 30y75p is the 75th percentile of the annual maximum temperatures recorded from 1992 to 2021.

Thermal Anomalies Analysis
For the evaluation of thermal anomalies, the median LST data of the four selected scenes were used as the reference surface temperature.To identify the areas with high LST values, the temperature difference (∆T) was calculated between the median LST value of each pixel and the average value of pixels ('pins') located in non-urban areas and 15-20 km away from the centroid of Padua.
The buffer size of 15-20 km used was defined according to the typical configuration of land cover/land use of the study area, which is the FUA.The ∆T was defined as the 'urban thermal anomaly' value.The main steps along with the tools used in QGIS are described below: 1.
Creation of a circular crown (15-20 km) 15 km distant from the FUA with an extension of 5 km.

2.
Extraction of the CLC layer using a 15-20 km circular crown level.

3.
From the extracted CLC database (15-20 km), identification of the points the maximum distance from the edge of the polygons classified as 'vegetation'.

4.
Selection of points at least 400 m away for the 15-20 km circular crown.

5.
Extraction of the value from the LST raster corresponding to the identified 'pins'.6.
Calculation of an average temperature value of the 'pins', which was then used to evaluate the thermal anomaly.

Relationship between Land Cover and Satellite Data
As previously mentioned, a raster database with a spatial resolution of 1 km was organised in order to overlay information from different sources (i.e., land cover and satellite data).The statistical analysis was carried out by overlying median LST and NDVI values (processed from four scenes from Sentinel-3), with the percentage of green areas extracted from the CLC database (resampled with a spatial resolution of 1 km).Using the GRASS tool 'r.regression.line' in QGIS, two correlations were performed, the first between the LST and percentage of green areas and the second one between the NDVI and percentage of green areas.

Heat-Related Urban Risk Assessment
Areas most at risk during HWs are those with high levels of thermal anomalies and a higher population density than other areas.Further, the most vulnerable people are the elderly (over 65) and children (younger than 5).Based on the literature [54][55][56][57][58][59], a GIS-based method was defined to identify urban and peri-urban areas in Padua most at risk.In particular, we focused on the elderly (the most vulnerable) by analysing the summer period (June-August).The HERI was calculated and mapped for the 31 municipalities located in the FUA by combining climate risk data-the thermal anomaly-and population data (i.e., population density and people over 65 years of age).

Exposure and Vulnerability Parameters
The exposure (E) and vulnerability (V) parameters were calculated for all municipalities within the FUA.The exposure parameter was assessed by computing the population density (inhabitants per km 2 ).The vulnerability (%) parameter was assessed by computing the percentage of the population aged 65 or over, the population most vulnerable to HWs.
The population data were collected from the statistical database of the Veneto Region (updated to 2020) at the municipal scale.Therefore, E and V values were calculated for each municipality.

Hazard Parameter
The hazard (H) parameter was obtained using the thermal anomaly values determined using the methodology presented in Section 2.3.The data included hourly air temperatures recorded by a weather station located near Padua, in Legnaro, and to four satellite images.The satellite images covered the following days: 1 August 2018, 27 June 2019, 31 July 2020 and 14 August 2021 (sensing hour UTC: 9.00-9.45a.m.).The images had a spatial resolution of 1 km 2 per pixel.Therefore, an average thermal anomaly value was calculated for each municipality.

Heat-Related Elderly Risk Index (HERI)
To quantify and map the HERI, the parameters E, V and H were normalised from 0 to 1 to enable comparison (E norm , V norm , H norm ).The normalisation procedure followed the 'max-min method' [60], which is described by Equation (4): where: y i is the normalised value in the dataset, which varies between 0 and 1 (E norm , V norm , H norm ); -x i is the value in the dataset; -min(x) is the minimum value in the dataset; -max(x) is the maximum value in the dataset.
Next, these values were combined into one index, the HERI, by assigning a weight to each parameter.Specifically, the parameters were combined by weighting E and V as 25%, respectively, and H as 50%: Finally, the HERI was normalised (HERI norm ) from 0 to 1 (always using the max-min method) and mapped at the municipal level.

Results
The results section is divided into five subsections.The first section presents the results of the UHI analysis, the identification of HW periods and the selection of satellite images.In Section 2, the indicators used to investigate the spatiotemporal distribution of the UHI phenomena are mapped.The third section describes the correlations between land cover and satellite data.The last two sections show the mapping of thermal anomalies and the HERI.

Urban Heat Island (UHI) Analysis
Based on a previous work [51], the HW periods were identified, and two indicators were calculated to investigate the UHI intensity.Subsequently, the satellite images were selected for the identified hot days based on the availability and quality of the Sentinel-3 images.
Table 2 lists the hottest days for each year (from 2016 to 2021), which also fit into an HW time window.Satellite images were processed for these selected days.Due to an acquisition issue in the years 2016 and 2017, this work analysed images from the following years: 2018, 2019, 2020 and 2021.For each year, a satellite image was selected and elaborated.
Since Sentinel-3 satellite images are acquired from two satellites, there is a greater availability of data than Landsat satellite images.The two in-orbit Sentinel-3 satellites enable a short revisit time of less than two days for OLCI and less than one day for SLSTR at the equator.

Mapping Urban Heat Waves and Heat Islands
The main data processed from the Sentinel-3 satellite images and from the CLC database were mapped for the urban and peri-urban areas of Padua, Italy, for summers of 2018-2021 (June-August).Using the QGIS 'resampling' tool, the median LST and NDVI values were calculated for the four selected satellite images.In Appendix A, the LST (median values of analysed scenes) is mapped with a nominal spatial resolution of 1 km 2 per pixel (Figure A1) and at a municipal level (Figure A2).

Normalised Difference Vegetation Index (NDVI) Calculation
The NDVI was mapped using a median NDVI value based on four satellite images (the same procedure used for the LST).From Figure 4, it can be seen that the lowest NDVI values (0.3-0.2), which correspond to areas with scarce or no vegetation, are located in the urban area of the FUA (see Figure A3, in Appendix A for a spatial resolution of 1 km 2 per pixel).

Urban Green Analysis: Estimation of Green Surfaces
Spatial analysis showed that higher NDVI values (0.71-0.61) occurred mainly in areas located in outer sectors of the FUA, in the municipalities of Teolo and Saccolongo (eastern sector) and Villafranca Padovana (northern sector).On the contrary, the lowest NDVI values (0.21-0.25) were mainly distributed within the city of Padua.
By estimating green surfaces on 1 km 2 area units, it is possible to geovisualise the interactions among land cover (Figure 5) and surface temperature (Figure 3), where high LST values and low NDVI values correspond to areas with a low percentage of green surfaces.

Correlations between Satellite Data and the Percentage of Green Surfaces
As described in Section 2.4, the statistical analysis was carried out in QGIS using the GRASS tool called 'r.regression.line'.Based on Equation (6), Table 3 shows the linear regression results between the median LST and NDVI values of the four satellite images selected and the percentage of green surfaces.These results show that there is a correlation between the LST (R is equal to −0.67) and the UHI effect and the presence of vegetation.As expected, the LST also depends on other factors that were not considered in this first investigation.However, this liner regression model can be used to evaluate how an increase in vegetation could mitigate the UHI effect by reducing the LST in the FUA of Padua.Correlating the percentage of green surfaces with the NDVI resulted in an R value of 0.76.

Mapping Urban Thermal Anomalies
Based on the analysis of thermal anomalies, the most critical municipality is Padua (ranked first with an average anomaly of 2.9 • C), followed by Ponte San Nicolò (2.5 • C anomaly), Noventa Padovana and Albignasego (2.4 • C anomaly) and Abano Terme (2.3 • C).
In the municipality of Padua, the most critical areas are the historic centre in the north part of the city, which is a densely built-up district with few green areas, and the industrial zone in the northeast of the municipality.These areas have the highest thermal anomaly values, equal to about 6-7 • C.Meanwhile, the lowest thermal anomaly values are in the peri-urban areas of Padua, where the building density is lower and the percentage of greenery is significantly higher than in the urban area.The areas least at risk are located in the municipality of Teolo, with minimum thermal anomaly values of −5.5 • C. In fact, Teolo is a municipality entirely included in the Regional Park of the 'Colli Euganei', dominated by forest in a gently hilled area.Urban thermal anomalies are shown in Figure 6.In Appendix A, urban thermal anomalies in the FUA are mapped with a spatial resolution of 1 km 2 per pixel (Figure A4) and at a municipal level (Figure A5).

Mapping the Heat-Related Elderly Risk Index (HERI)
The results of this analysis show that the most vulnerable municipality during HWs is Padua, while the least vulnerable is Polverara.
In general, the municipalities with the highest levels of population density (defined as exposure, E) are Padua (2217 inhabitants per km 2 ) and Noventa Padovana (1537 inhabitant per km 2 ).Meanwhile, the municipalities with the lowest E values are Bovolenta and Terrassa Padovana, with a population density of less than 200 inhabitants per km 2 .
The highest thermal anomaly values (defined as hazard, H) during the summer period between 2018 and 2021 were observed in Padua (2.91 • C on average), followed by Ponte San Nicolò (2.54 • C on average), while the lowest H values were observed in Teolo (−0.96 • C on average) and Poverara (−0.3 • C on average).
Table 4 presents the non-normalised values of E, V, H and HERI for each municipality.The municipalities are sorted according to HERI values, from lowest to highest. Figure 7 shows the HERI values at the municipal level according to five risk levels: very low (HERI norm ≤ 0.2), low (0.2 < HERI norm ≤ 0.4), moderate (0.4 < HERI norm ≤ 0.6), high (0.6 < HERI norm ≤ 0.8) and very high (HERI norm > 0.8).The very high-risk levels are localised in the urban area of Padua and in three municipalities near the city (Abano Terme, Noventa Padovana and Ponte San Nicolò).These municipalities represent 19.8% of the FUA area, equal to 121.5 km 2 .On the contrary, 16 municipalities have low and very low HERI levels, representing 48.1% (395.4 km 2 ) of the total area.In Appendix A, the three parameters used to identify the most vulnerable areas (i.e., E, V and H) are mapped in Figures A6-A8.

Discussion
Several studies have shown that there are different methodologies and techniques for investigating and mapping summer thermal anomalies in urban areas.GIS tools are fundamental for carrying out this type of analysis, and there are several variables that can be processed and mapped.For example, Guerri et al. [61] assessed different indicators able to describe the land cover (e.g., LST, NDVI and albedo), the urban morphology and the distribution of the population.The main input data are satellite images, and based on the analysis scale, the type of resolution and the availability of the data, Landsat or Sentinel images can be used [62,63].
In this work, Sentinel-3 satellite images were processed and the main indicators elaborated to investigate the urban thermal anomalies are LST, NDVI and HERI.The results of this study show that the areas most at risk are those with LST values higher than 40 • C, NDVI lower than 0.4 and HERI higher than 0.6-0.65.In particular, HERI values higher than 0.8 identify high-risk areas, while values between 0.6 and 0.8 indicate areas at risk.
The UHI phenomenon has important impacts on public health.Europe is highly vulnerable to HWs [64].UHI mitigation measures should be adopted to reduce the mortality and vulnerability of older people as well as other weak stakeholders (i.e., poor people, homeless).Mitigation measures include increasing urban green spaces, installing cool roofs (e.g., reflective roofs) and green walls, updating building regulations to promote the use of passive cooling strategies and introducing financing instruments (i.e., incentives) to support the cost of the interventions [65].Such measures should be integrated into climate-resilient development plans that consider both the physical impacts of UHI on urban territory and the spatial dimensions of the most vulnerable people.In fact, two of the main issues related to adaptation strategies are their effectiveness and spatial inequality.Currently, there is no mechanism for monitoring and reporting adaptation strategies, which would allow for the supervision and evaluation of the adaptation process and support the involvement of weak actors in urban planning strategies [66,67].Regarding spatial inequality, the mitigation measures might be localised in a selected area of the urban territory to the detriment of the most socially disadvantaged areas.In fact, social inequality affects the localisation of adaptation and mitigation strategies.To overcome this inequality, cities should promote inclusive sustainable development [7,68].For example, encouraging the spatial equity of urban green spaces would be an adoptable solution to mitigate the UHI effect [69].Of the various aspects that influence the UHI phenomenon, the urban form plays an essential role.Thus, as suggested by Mohtat and Khirfan [70], it is necessary to integrate climate justice into the adaptation of the urban form.The density and size of cities affect the UHI intensity [71].More specifically, the urban parameter that has the greatest influence on LST values is the green space ratio, followed by the building coverage ratio, the building ratio and the height-to-width ratio [72,73].

Conclusions
To mitigate UHI and address the problem of climate change, it is necessary to promote sustainable and resilient development.By identifying the most at-risk areas, targeted mitigation measures could be designed to improve the outdoor thermal conditions and liveability of urban spaces.

Main Research Findings
The present work involved mapping thermal anomalies in the urban and peri-urban areas of Padua, Italy, (FUA geographical scale) using an integrated GIS and remote-sensing technique based on the processing of satellite images.Due to the scale of the analysis and the availability of data, Sentinel-3 images were used to assess the UHI and map the LST, NDVI, percentage of green surfaces, urban thermal anomalies and the HERI at different scales and raster resolutions.
The results showed that the most critical areas with higher temperatures are in the municipality of Padua, with maximum thermal anomaly values of approximately 6-7 • C.An analysis on urban heat-related health risks was also carried out.The HERI was mapped at the municipal level.Hazardous HERI levels (high and very high risk) were localised in Padua, Abano Terme, Noventa Padovana and Ponte San Nicolò.
Future studies could define adaptation strategies to be applied in the most vulnerable areas and assess their impacts on the territory in terms of outdoor thermal comfort and energy demand for cooling.

Limitations and Future Direction of the Research
The use of satellite data such as Sentinel-3 provides an effective open-access and low-cost monitoring system to assess and to monitor extreme heat wave events and UHI intensity on urban territory at multiple scales over time.In fact, as reported in the literature, air temperature data with a suitable spatial distribution in urban areas of WMO certified ground weather stations are often not available.In fact, ground weather stations provide air data temperature only in a specific point of the territory and they might be costly if not supported by public institutions.Therefore, even with some limitations, in the case of a lack of spatially consistent weather stations, the use of LST derived from Sentinel-3 data provides, during extreme heat wave events, a synoptic assessment of UHI intensity in wide portions of the territory (local and supra-local geographical scales).Moreover, thanks to the high-frequency re-visit time of Sentinel-3, further research should focus on the spatial and temporal variability of thermal anomalies on the urban territory to give precise and practical indications for climate resilient urban planning and development.

Figure 2 .
Figure 2. Flowchart of the GIS-based methodology.

3. 2 . 1 .
Mapping Land Surface Temperature (LST) Spatial analysis of LSTs showed that maximum LST values (44.5-45 • C) are localised within the municipality of Padua (Figure 3).The areas most affected by UHI are mainly located in the industrial zone (eastern sector) and in the central-northern neighbourhoods of Padua.Here, LST values ranges from about 41 to 43 • C. High LST values are also localised in the NW sector, outside the city, in the municipalities of Rubano and Mestrino.

Figure 3 .
Figure 3. Map of the land surface temperature in the Functional Urban Area of Padua (NE Italy) based on four Sentinel-3 scenes (100 m raster resolution output).

Figure 4 .
Figure 4. NDVI map (median value) geovisualising green areas in the Functional Urban Area of Padua (100 m raster resolution output).

Figure 5 .
Figure 5. Map of the percentage of green surfaces with a 1 km 2 grid size in the Functional Urban Area of Padua in Italy using the Corinne Land Cover (CLC) in 2018.

Figure 6 .
Figure 6.Map of urban thermal anomalies in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (100 m raster resolution output).

Figure 7 .
Figure 7. Map of the heat-related elderly risk index (HERI) in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (normalised values from 0 to 1).

Figure A2 .
Figure A2.Map of LSTs at the municipal level in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (June-August).

Figure A3 .
Figure A3.NDVI map (median value) geovisualising green areas in the Functional Urban Area of Padua (1 km raster resolution output).

Figure A4 .
Figure A4.Map of urban thermal anomalies in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (1 km raster resolution output).

Figure A5 .
Figure A5.Map of urban thermal anomalies at the municipal level in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (June-August).

Figure A6 .
Figure A6.Map of the exposure (E) parameter in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (normalised values from 0 to 1).

Figure A7 .
Figure A7.Map of the vulnerability (V) parameter in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (normalised values from 0 to 1).

Figure A8 .
Figure A8.Map of the hazard (H) parameter in the Functional Urban Area of Padua in Italy during the summers of 2018-2021 (normalised values from 0 to 1).

Table 1 .
Description of spatial input data.

Table 2 .
Identification of Sentinel-3 satellite images according to HW time windows.
* Satellite images not applicable due to acquisition issue.

Table 3 .
Linear regression results between satellite data (LST and NDVI) and the percentage of green surfaces.
1 R is the correlation between the predictor variable and the response variable.

Table 4 .
Non-normalised values of E, V, H and HERI in the urban and per-urban areas of Padua in Italy during the summers of 2018-2021.