Response of Urban Heat Stress to Heat Waves in Athens (1960–2017)

: The increasing frequency, intensity and duration of heat waves seem to follow the observed global warming in recent decades. Vulnerability to heat waves is expected to increase in urban environments mainly due to population density and the e ﬀ ect of the urban heat island that make cities hotter than surrounding non-urban areas. The present study focuses on a vulnerable area of the eastern Mediterranean, already characterized as a ‘hot spot’ with respect to heat-related risk and investigates the change in heat stress levels during heat wave compared to non-heat wave conditions as well as the way that heat stress levels respond to heat waves in urban, compared to non-urban, environments. The adoption of a metric accounting for both the intensity and duration of the hot event yielded a total of 46 heat wave episodes over a nearly 60-year period, but with very rare occurrence until the late 1990s and a profound increased frequency thereafter. The results reveal a di ﬀ erence of at least one thermal stress category between heat wave and non-heat wave periods, which is apparent across the entire range of the thermal stress distribution. The analysis demonstrates a robust intensiﬁcation of nighttime heat stress conditions in urban, compared to non-urban, sites during severe heat waves. Nevertheless, severe heat waves almost equalize heat stress conditions between urban and non-urban sites during midday.


Introduction
Unlike other climate hazards, heat waves (hereafter HWs) have only lately been acknowledged as a major threat to environment and society [1], although they are related to more deaths than other extreme events [2]. The increasing frequency, intensity and duration of HWs in different world regions seem to follow the observed global warming in recent decades [3,4]. The impact of HWs on human health has been investigated in terms of mortality, morbidity and hospital admissions as well as in terms of heat stress and human body's energy balance [5][6][7][8][9][10][11][12]. The populations in low and middle income countries are more affected from HWs as compared to high income countries and as global warming progresses, they will still be more affected [13].
A period of consecutive days with exceptionally hot weather which deviates from the local climatic conditions is usually considered as a HW. Nevertheless, there is no standard and universal definition for HWs. Several metric systems (indices) and approaches have been proposed in order to detect and investigate these abnormal heat events based on either absolute or relative (percentile) thresholds. Most of the approaches try to portray the frequency, the amplitude and the persistence of HWs mainly based on expressions of air temperature, such as the maximum, the minimum and the daily mean (e.g., [4,14,15]). Despite the fact that HWs are meteorological events associated with large scale atmospheric anomalies and persistent anticyclonic conditions, they should also be assessed in terms Given the fact that the present study aimed to identify HWs in a large urban agglomeration and assess the heat stress risk on its population during these periods, an appropriate station with long, continuous and high quality records was required. National Observatory of Athens (NOA) met these preconditions and the station has received a formal recognition by the World Meteorological Organization (WMO) for its centennial operation. NOA is located on a hill (107 m a.s.l.) at the Athens urban core. The air temperature time series have been quality controlled for their homogeneity [39] whilst the long-term surface solar radiation time series have been quality controlled for their accuracy [40].
The daily maximum air temperature (T max ) time series were employed for a long 58-year period from 1960 to 2017. Moreover, the hourly time series of air temperature, relative humidity, wind speed and global solar radiation were also used for the aforementioned time period. The identification of abnormal heat events was focused on the warm period of the year and hence, the used data were restricted to the summer months of June, July and August. No missing values were found for the T max time series, whereas less than 0.6% in total were detected for the hourly time series of air temperature, relative humidity, wind speed and global solar radiation.
In order to investigate the spatial differences of heat stress conditions between an urban and a suburban area during exceptional abnormal heat events, the suburban station of Tatoi was also employed. Tatoi station (237 m a.s.l.) is located at the north boundaries of the Athens urban complex between the mountains of Parnitha and Hymettus and is 16 km distant from the NOA station. The 3-hourly time series of air temperature, relative humidity, wind speed and cloud cover were also employed for the 1985-2010 period.

Heat Wave Identification Process
The identification of HWs was attained employing a day-count indicator based on percentile thresholds in a span of at least three consecutive days. In particular, the period where the daily T max exceeded the 95th percentile of the reference period of 30 years (1971-2000) for more than three consecutive days was characterized as a 'heat wave' period (HW period). It is known that literature lacks a universal definition of HWs (e.g., [16]). The use of upper percentiles in T max and the choice of the least duration of three consecutive days in the adopted definition is consistent (or more strict) with definitions used in other relevant studies [41][42][43][44].
The 95th percentile of the reference period was computed from five-day windows centered on each calendar day. The selection of a five-day window approach produced a sample size of 150 values (30 years × 5 days) on each calendar day. Such an approach expressed the anomalies relative to the local climatic conditions and created a relatively smooth cycle of percentiles thresholds for each calendar day of summer based on the reference period   [14]. According to the WMO [14], the selection of another reference period had only a small impact on the results over time. This procedure could be applied to all seasons creating an annual cycle of percentile thresholds but other periods could also be selected. As already mentioned, the analysis was seasonally restricted and hence, the results were focused on the summer season where the impact of HWs on thermal comfort of urban population was maximized. In order for the method to be applicable to the first days of June and the last days of August, the T max time series of the last days of May and the first days of September were also employed.
In order to assess the intensity and severity of HWs, two additional metrics were used. 'Cumulative exceedance' (CumExc) was defined as the sum of the differences between the T max and the 95th T max percentile for all days of the HW episode. 'Absolute exceedance' (AbsExc) was defined as the maximum daily T max exceedance above the 95th T max percentile during the HW episode.
where i indicates a calendar day of the HW and T max,P95 denotes the 95th T max percentile of the day.

Heat Stress Assessment
The assessment of heat stress was performed by employing one simple 'heat-humidity' index, the so-called Humidex (HD), and one thermo-physiologically based energy balance index, the so-called Universal Thermal Climate Index (UTCI). HD is an empirical index based on an algebraic/statistical model, omitting significant energy fluxes, while UTCI is an energy balance stress index [20][21][22].
HD follows a simplistic approach to assess the degree of comfort taking into consideration the combined effect of air temperature and humidity on physical distress [45]. However, HD index does not take into account the parameters of radiation and wind speed, omitting also behavioral factors such as the metabolic rate and the clothing insulation as well as the energy balance equation between the human body and the thermal environment [22]. This index has been exclusively designed for warm weather conditions and thus, its result is valid when the ambient temperature is higher than 21.0 • C ( Table 1). The HD calculations were performed following the formula: where T is the air temperature ( • C) and e is the vapor pressure (hPa). According to the guide for summer weather hazards released by the Environment Canada Service, an HD value higher than 40 can be considered as an extremely high HD reading, where people should curtail all unnecessary outdoor activities [46]. In addition, certain types of activity should be reduced when the HD value is in the mid to high 30s. It is also worth noting that the assessment scale of HD has been calibrated for people living in temperate regions of Canada. Table 1. Degree of comfort according to Humidex (HD) [46].

>45
Dangerous-heat stroke possible [40][41][42][43][44][45] Great discomfort-avoid exertion [30][31][32][33][34][35][36][37][38][39] Some discomfort ≤29 No discomfort The UTCI is derived from an energy balance model that takes into account both the mechanisms of heat exchange and the thermo-physiological processes on the human body. Specifically, UTCI is based on an advanced multi-node model of thermoregulation [47] which is coupled with an adaptive clothing model [48]. UTCI is expressed as equivalent temperature of a reference environment and its applicability covers a wide range of environmental conditions. Table 2 presents the thermal stress categories according to UTCI [49]. Unlike the HD index, whose input variables are only the air temperature and humidity, the meteorological input variables of UTCI include air temperature, air humidity, wind speed and mean radiant temperature. The adaptive clothing model [48] integrated in UTCI-Fiala model [47] determines the clothing insulation, while the metabolic rate is set at 135 W·m −2 (2.3 met) and the walking speed at 1.1 m·s −1 [50]. Hence, UTCI describes the thermal environment of human beings in a more comprehensive way and thus, is a more robust and advanced index compared to HD [20,21]. According to Alfano et al. [22], the HD index shows significant shortcomings and thus, it is not suitable for assessing heat stress both outdoors and indoors. Nevertheless, the HD index is still used by researchers or National Weather Services for communication purposes regarding heat-related hazards, overlooking its shortcomings due to its simplicity of calculations.
The UTCI calculations were performed based on a Fortran subroutine retrieved from the official site of UTCI project (http://www.utci.org/). The simulations of radiant temperature, which is one of the input variables of UTCI, were carried out employing the RayMan model that simulates long-wave and short-wave radiation fluxes [51,52].
Following the definition of 'Urban Heat Island Intensity' (UHII) which is quantified by the difference of air temperature between urban and non-urban areas [23], we further proposed the indicator of 'Urban Heat Stress Intensity' (UHSI) which was quantified as the difference of heat stress between the urban station of NOA and the suburban station of Tatoi. Hence, the UHSI is given by the following equation: where the 'Metric' indicates the bioclimatic index (HD or UTCI) used for the assessment of heat stress intensity between urban and suburban areas.

Heat Wave Identification
In order to identify the HWs, the 95th T max percentile-based thresholds on each calendar day of summer were produced for the reference period (1971-2000) ( Figure 1). The percentile thresholds varied during summer, with the values ranging between 31.8 • C and 39.0 • C. The lowest values of thresholds were observed in early June and the highest values were noted during July as they were associated with the local climatic conditions. Atmosphere 2019, 10, x FOR PEER REVIEW 5 of 17 radiant temperature. The adaptive clothing model [48] integrated in UTCI-Fiala model [47] determines the clothing insulation, while the metabolic rate is set at 135 W·m −2 (2.3 met) and the walking speed at 1.1 m·s −1 [50]. Hence, UTCI describes the thermal environment of human beings in a more comprehensive way and thus, is a more robust and advanced index compared to HD [20,21]. According to Alfano et al. [22], the HD index shows significant shortcomings and thus, it is not suitable for assessing heat stress both outdoors and indoors. Nevertheless, the HD index is still used by researchers or National Weather Services for communication purposes regarding heat-related hazards, overlooking its shortcomings due to its simplicity of calculations. The UTCI calculations were performed based on a Fortran subroutine retrieved from the official site of UTCI project (http://www.utci.org/). The simulations of radiant temperature, which is one of the input variables of UTCI, were carried out employing the RayMan model that simulates long-wave and short-wave radiation fluxes [51,52].
Following the definition of 'Urban Heat Island Intensity' (UHII) which is quantified by the difference of air temperature between urban and non-urban areas [23], we further proposed the indicator of 'Urban Heat Stress Intensity' (UHSI) which was quantified as the difference of heat stress between the urban station of NOA and the suburban station of Tatoi. Hence, the UHSI is given by the following equation: where the 'Metric' indicates the bioclimatic index (HD or UTCI) used for the assessment of heat stress intensity between urban and suburban areas.

Heat Wave Identification
In order to identify the HWs, the 95th Tmax percentile-based thresholds on each calendar day of summer were produced for the reference period (1971-2000) ( Figure 1). The percentile thresholds varied during summer, with the values ranging between 31.8 °C and 39.0 °C. The lowest values of thresholds were observed in early June and the highest values were noted during July as they were associated with the local climatic conditions.   Figure 2 shows the frequency of HWs per year for the study period 1960-2017. The methodological approach used identified 46 HWs in total, during the long 58-year period. Only three HWs were detected until the mid-1980s, while the HWs frequency increased markedly afterwards. The results revealed an average number of 0.8 HWs y −1 for the entire study period. However, the average number of HWs was noticeably larger in the 2000s and 2010s, at the level of 1.9 HWs y −1 and 1.4 HWs y −1 , respectively. The highest number of HWs was observed in 2007 and 2012, both with five episodes, and the year 2000 followed with four episodes. The frequency of HWs per month indicated that 21.7% (10 events) occurred in June, 28.3% (13 events) occurred in July and 43.5% (20 events) occurred in August. In addition, 6.5% of HWs (three events) took place between two months (June/July and July/August).
Atmosphere 2019, 10, x FOR PEER REVIEW 6 of 17 Figure 2 shows the frequency of HWs per year for the study period 1960-2017. The methodological approach used identified 46 HWs in total, during the long 58-year period. Only three HWs were detected until the mid-1980s, while the HWs frequency increased markedly afterwards. The results revealed an average number of 0.8 HWs y −1 for the entire study period. However, the average number of HWs was noticeably larger in the 2000s and 2010s, at the level of 1.9 HWs y −1 and 1.4 HWs y −1 , respectively. The highest number of HWs was observed in 2007 and 2012, both with five episodes, and the year 2000 followed with four episodes. The frequency of HWs per month indicated that 21.7% (10 events) occurred in June, 28.3% (13 events) occurred in July and 43.5% (20 events) occurred in August. In addition, 6.5% of HWs (three events) took place between two months (June/July and July/August).  Figure 3 depicts the CumExc and AbsExc per HW as well as the duration and the year of appearance for each HW, organized in descending order of CumExc. The duration of the 46 HWs ranged between 3 to 11 days while 18 HW events spanned a period of 3 days. Only one HW was detected with an 11 day duration followed by five HWs with an 8 day duration. The CumExc ranged from 0.3 °C to 26.7 °C and the AbsExc varied between 0.2 °C and 7.8 °C. The HW with the highest CumExc displayed both the longest duration and the highest AbsExc. The second ranked HW based on CumExc was also second in terms of duration and third in terms of AbsExc. Nevertheless, a visual inspection of the results revealed significant differences in CumExc and AbsExc even in HWs with the same duration. It is also noteworthy that the first 15 HWs according to the ranking based on CumExc all originated after 1999, with the exception of 2 HW events.  Figure 3 depicts the CumExc and AbsExc per HW as well as the duration and the year of appearance for each HW, organized in descending order of CumExc. The duration of the 46 HWs ranged between 3 to 11 days while 18 HW events spanned a period of 3 days. Only one HW was detected with an 11 day duration followed by five HWs with an 8 day duration. The CumExc ranged from 0.3 • C to 26.7 • C and the AbsExc varied between 0.2 • C and 7.8 • C. The HW with the highest CumExc displayed both the longest duration and the highest AbsExc. The second ranked HW based on CumExc was also second in terms of duration and third in terms of AbsExc. Nevertheless, a visual inspection of the results revealed significant differences in CumExc and AbsExc even in HWs with the same duration. It is also noteworthy that the first 15 HWs according to the ranking based on CumExc all originated after 1999, with the exception of 2 HW events.
The first two ranked HWs based on CumExc are two exceptional, and at the same time representative, abnormal heat events of Greece. The first one took place in June 2007 and the second in July 1987 with significant implications for the population, especially for the inhabitants of densely urban areas like Athens [34,[53][54][55][56][57][58].

Changes in Heat Stress between HW and Non-HW Periods
It is interesting to assess and quantify the magnitude of the heat stress change between HW and non-HW periods. To this end, the cumulative distribution function (CDF) curves of the bioclimatic indices (HD and UTCI) were produced during HW and non-HW periods ( Figure 4). CDF curves are useful to find probabilities above or below a certain value of the index. The period of HWs refers exclusively to the 46 events of the 1960-2017 period, while the non-HW period refers to the reference period  excluding the detected HWs of that period. The results were based on the hourly time series of the selected indices over these periods. Hence, the CDF curves were the result of applying a number of bins across the range of HD and UTCI hourly time series and estimating the corresponding frequencies. The first two ranked HWs based on CumExc are two exceptional, and at the same time representative, abnormal heat events of Greece. The first one took place in June 2007 and the second in July 1987 with significant implications for the population, especially for the inhabitants of densely urban areas like Athens [34,[53][54][55][56][57][58].

Changes in Heat Stress between HW and Non-HW Periods
It is interesting to assess and quantify the magnitude of the heat stress change between HW and non-HW periods. To this end, the cumulative distribution function (CDF) curves of the bioclimatic indices (HD and UTCI) were produced during HW and non-HW periods ( Figure 4). CDF curves are useful to find probabilities above or below a certain value of the index. The period of HWs refers exclusively to the 46 events of the 1960-2017 period, while the non-HW period refers to the reference period  excluding the detected HWs of that period. The results were based on the hourly time series of the selected indices over these periods. Hence, the CDF curves were the result of applying a number of bins across the range of HD and UTCI hourly time series and estimating the corresponding frequencies.
According to Figure 4, the CDF curves of both indices are shifted noticeably towards higher values of the index during the HW period. This shift is observed not only in the upper but also in the lower levels of the index, which indicates that the heat stress conditions during HWs are significantly worse at nighttime as well, when lower values of the index were observed. The difference between HW and non-HW periods was one or even two stress categories of the index scale, with the exception of percentile ranks at the ends of the distribution. Therefore, significant percentile changes between non-HW and HW periods could be observed across a wide range of UTCI and HD values. In particular, the HD values from 25.5 to 42.5 showed changes of at least 10.0 percentage points between non-HW and HW periods, and this difference reached up to 47.4 percentage points at the 32.5 scale point (some discomfort) (Figure 4a). The results also reveal that the frequency of cases with HD ≥ 40, which indicates at least 'great discomfort-avoid exertion' (Table 1), rose up to 23.0% during HWs compared to the frequency of 1.4% which corresponded to non-HW conditions. Regarding the 'dangerous-heat stroke possible' category (HD > 45), the frequency reached up to 3.6% in HWs, while it was negligible (0.03%) in non-HW conditions.  Table 1 and Table 2.
For the UTCI, a wider range of values was observed with significant percentile changes between non-HW and HW periods compared to the corresponding range of HD values (Figure 4b). Therefore, UTCI differed by at least 10.0 percentage points in the value range between 17.7 °C and 41.4 °C and the difference reached up to 37.1 percentage points at the first observed peak (23.7 °C) and up to 25.2 percentage points at the second observed peak (35.7 °C). The frequency of cases with UTCI > 38 °C, which indicates at least very strong heat stress (Table 2), reached up to 25.2% in HW According to Figure 4, the CDF curves of both indices are shifted noticeably towards higher values of the index during the HW period. This shift is observed not only in the upper but also in the lower levels of the index, which indicates that the heat stress conditions during HWs are significantly worse at nighttime as well, when lower values of the index were observed. The difference between HW and non-HW periods was one or even two stress categories of the index scale, with the exception of percentile ranks at the ends of the distribution. Therefore, significant percentile changes between non-HW and HW periods could be observed across a wide range of UTCI and HD values. In particular, the HD values from 25.5 to 42.5 showed changes of at least 10.0 percentage points between non-HW and HW periods, and this difference reached up to 47.4 percentage points at the 32.5 scale point (some discomfort) (Figure 4a). The results also reveal that the frequency of cases with HD ≥ 40, which indicates at least 'great discomfort-avoid exertion' (Table 1), rose up to 23.0% during HWs compared to the frequency of 1.4% which corresponded to non-HW conditions. Regarding the 'dangerous-heat stroke possible' category (HD > 45), the frequency reached up to 3.6% in HWs, while it was negligible (0.03%) in non-HW conditions.
For the UTCI, a wider range of values was observed with significant percentile changes between non-HW and HW periods compared to the corresponding range of HD values (Figure 4b). Therefore, UTCI differed by at least 10.0 percentage points in the value range between 17.7 • C and 41.4 • C and the difference reached up to 37.1 percentage points at the first observed peak (23.7 • C) and up to 25.2 percentage points at the second observed peak (35.7 • C). The frequency of cases with UTCI > 38 • C, which indicates at least very strong heat stress (Table 2), reached up to 25.2% in HW period, while the corresponding frequency under non-HW conditions was only 2.6%, indicating a tenfold increase during HWs. In addition, the frequency of cases with UTCI < 26 • C, which suggests no thermal stress, rose up to 57.3% under non-HW conditions compared to the frequency of 27.3% corresponding to the HW period. Figure 5 presents the mean diurnal variation of thermal stress levels as reflected in UTCI and HD indices during the HW and non-HW periods. The HD results show that 'great discomfort-avoid exertion' (HD ≥ 40) was observed between 11:00 and 17:30 local standard time (LST) during the HW period (Figure 5a). On the contrary, there was an absence of 'great discomfort-avoid exertion' under non-HW conditions and 'some discomfort' (30 ≤ HD < 40) was noted only from 09:00 to 20:00 LST. Although 'some discomfort' was observed during daytime under non-HW conditions, the nighttime thermal conditions always indicated 'no discomfort'. This is particularly important as it suggests that during nighttime the human body can recover from daytime heat stress. As far as the UTCI results are concerned, very strong heat stress was observed between 10:00 and 17:00 LST while strong heat stress was noted from 08:00 to 11:00 LST and from 17:00 to 19:00 LST during the HWs period (Figure 5b). On the contrary, the very strong heat stress category was not recorded under non-HW conditions and only the category of strong heat stress was noted between 10:00 and 16:00 LST.
The diurnal range of both UTCI and HD was higher during the HW period than under non-HW conditions ( Figure 5). Thus, the diurnal range of UTCI (HD) was 17.6 • C (10.3) and 15.2 • C (7.1) during the HW and non-HW periods, respectively. It is also noted that the maximum value of mean heat stress (observed during midday hours) was about 7.4 • C and 7.3 scale points higher during HW periods than during non-HW periods for the UTCI and HD, respectively. The corresponding absolute increase for the minimum value of mean thermal stress (observed during nighttime) was approximately 5.0 • C for the UTCI and 4.1 scale points for the HD in HW periods compared to non-HW periods. Figure 6 presents the total exposure time in hours per HW when people experience heat stress, organized in descending order of CumExc. As already presented in Figure 3, the CumExc refers to the sum of differences between the T max and the 95th T max percentile threshold for all days of each HW. Based on HD results, the total exposure time when the population experienced at least 'some discomfort' (HD ≥ 30) ranged between 40 and 260 h for the detected HWs (Figure 6a). In addition, there were five HWs with no exceedance of the next upper threshold of HD which indicates at least 'great discomfort-avoid exertion' (HD ≥ 40) and thus, the exposure time for this threshold ranged from 0 to 98 h. The UTCI results show that the exposure time with at least strong and very strong heat stress varied from 22 to 128 hours and from 1 to 77 h, respectively (Figure 6b).  Table 1 and Table 2.
The diurnal range of both UTCI and HD was higher during the HW period than under non-HW conditions ( Figure 5). Thus, the diurnal range of UTCI (HD) was 17.6 °C (10.3) and 15.2 °C (7.1) during the HW and non-HW periods, respectively. It is also noted that the maximum value of mean heat stress (observed during midday hours) was about 7.4 °C and 7.3 scale points higher during HW periods than during non-HW periods for the UTCI and HD, respectively. The corresponding absolute increase for the minimum value of mean thermal stress (observed during nighttime) was approximately 5.0 °C for the UTCI and 4.1 scale points for the HD in HW periods compared to non-HW periods. Figure 6 presents the total exposure time in hours per HW when people experience heat stress, organized in descending order of CumExc. As already presented in Figure 3, the CumExc refers to the sum of differences between the Tmax and the 95th Tmax percentile threshold for all days of each HW. Based on HD results, the total exposure time when the population experienced at least 'some discomfort' (HD ≥ 30) ranged between 40 and 260 hours for the detected HWs (Figure 6a). In addition, there were five HWs with no exceedance of the next upper threshold of HD which indicates at least 'great discomfort-avoid exertion' (HD ≥ 40) and thus, the exposure time for this threshold ranged from 0 to 98 hours. The UTCI results show that the exposure time with at least strong and very strong heat stress varied from 22 to 128 hours and from 1 to 77 hours, respectively (Figure 6b).  The general fluctuating pattern of exposure time under heat stress followed to some extent the pattern of CumExc. As the CumExc decreased, the exposure time under heat stress declined as well, however, this was observed for almost the top half of the ranked HWs. Subsequently, the exposure time showed no clear downward trend. The CumExc depended to some extent on the duration of the HW which in turn affected the potential exposure time under heat stress, however, this was not always the case even in HWs with the same duration. Hence, there were HW events with the same duration and significant difference in CumExc but with similar number of hours under heat stress. For example, the HW event that ranked 13th based on CumExc (10.8 °C) had a duration of 4 days and commenced at the end of July 2000 (Figure 6b). During this HW, the number of hours exceeding the analyzed thresholds was 47 hours for UTCI > 32 °C and 33 hours for UTCI > 38 °C. On the contrary, the HW event that ranked 32nd and took place in early August 2014, showed considerably lower CumExc (3.6 °C) although it lasted longer (5 days), while the number of hours exceeding the analyzed thresholds was similar compared to the number of hours of the HW event that ranked 13th. This is attributed to the fact that the ranking of HWs was based on Tmax exceedances, however, heat stress does not only depend on air temperature as other atmospheric parameters (i.e., relative humidity, wind speed, radiation) play critical roles as well. This implies that even less severe HWs (in terms of intensity and duration) may pose severe heat stress to the population. The general fluctuating pattern of exposure time under heat stress followed to some extent the pattern of CumExc. As the CumExc decreased, the exposure time under heat stress declined as well, however, this was observed for almost the top half of the ranked HWs. Subsequently, the exposure time showed no clear downward trend. The CumExc depended to some extent on the duration of the HW which in turn affected the potential exposure time under heat stress, however, this was not always the case even in HWs with the same duration. Hence, there were HW events with the same duration and significant difference in CumExc but with similar number of hours under heat stress. For example, the HW event that ranked 13th based on CumExc (10.8 • C) had a duration of 4 days and commenced at the end of July 2000 (Figure 6b). During this HW, the number of hours exceeding the analyzed thresholds was 47 h for UTCI > 32 • C and 33 h for UTCI > 38 • C. On the contrary, the HW event that ranked 32nd and took place in early August 2014, showed considerably lower CumExc (3.6 • C) although it lasted longer (5 days), while the number of hours exceeding the analyzed thresholds was similar compared to the number of hours of the HW event that ranked 13th. This is attributed to the fact that the ranking of HWs was based on T max exceedances, however, heat stress does not only depend on air temperature as other atmospheric parameters (i.e., relative humidity, wind speed, radiation) play critical roles as well. This implies that even less severe HWs (in terms of intensity and duration) may pose severe heat stress to the population. Figure 7 presents the diurnal variation of UHSI as defined in Section 2.3 based on HD and UTCI indices for an exceptional HW and a non-HW period. The HW of June 2007 was selected as the exceptional HW, since it showed the greatest intensity and severity based on both the T max (Figure 3) and the frequency of hours under heat stress among the detected HWs ( Figure 6). The non-HW period was represented from the same calendar days that correspond to the HW of June 2007, covering the 1985-2010 period (following the data availability of the suburban station) and excluding any HWs that fall in these days. Such an approach expresses the local climatic conditions for the period when the exceptional HW of 2007 took place, excluding the abnormal heat events of this period. A positive UHSI indicates higher heat stress at the urban area compared to the suburban one. Striking differences of UHSI features between strong HW and non-HW conditions were observed depending on the time of the day (Figure 7). Based on HD results, the magnitude of UHSI held a positive sign of the order of 2.0 scale points with weak diurnal variation during non-HW conditions, suggesting higher heat stress conditions at the urban station throughout the 24-h period (Figure 7a). On the contrary, the diurnal variation of UHSI was markedly modified during the HW episode, showing a strong amplification during nighttime compared to non-HW conditions. The magnitude of UHSI exceeded 7.0 scale points during nighttime, being up to 5.0 scale points greater compared to non-HW conditions. However, UHSI declined significantly during the midday hours (11:00-16:00 LST) and at 14:00 LST, it marginally held a negative sign. This pattern indicates that intense HWs exacerbated urban heat stress during nighttime, while during midday both urban and suburban sites experienced comparable heat stress conditions. It is also worth mentioning that the diurnal pattern of UHSI during HWs resembled to some extent the typical UHII pattern with maxima during nighttime and minima during daytime.

Urban Heat Stress Intensity between HW and Non-HW Conditions
According to UTCI, the magnitude of UHSI revealed strong diurnal variation in non-HW conditions, with positive and higher values up to 4.2 °C occurring during late evening and nighttime hours, suggesting in general higher heat stress at the urban site (Figure 7b). UHSI responded to strong HWs by amplification during nighttime, when its magnitude reached 6.4 °C (being higher by at least 2.0 °C compared to non-HW conditions) and decreased during midday (11:00-16:00 LST). The mean diurnal pattern of the UHSI during the strong HW presented a kind of U-shaped curve, though more flattened than the corresponding steeper curve based on HD. This difference is probably attributed to the impact of insolation and wind on thermal stress, since these variables were incorporated in the UTCI calculations but they were omitted in the HD calculations. Striking differences of UHSI features between strong HW and non-HW conditions were observed depending on the time of the day (Figure 7). Based on HD results, the magnitude of UHSI held a positive sign of the order of 2.0 scale points with weak diurnal variation during non-HW conditions, suggesting higher heat stress conditions at the urban station throughout the 24-h period (Figure 7a). On the contrary, the diurnal variation of UHSI was markedly modified during the HW episode, showing a strong amplification during nighttime compared to non-HW conditions. The magnitude of UHSI exceeded 7.0 scale points during nighttime, being up to 5.0 scale points greater compared to non-HW conditions. However, UHSI declined significantly during the midday hours (11:00-16:00 LST) and at 14:00 LST, it marginally held a negative sign. This pattern indicates that intense HWs exacerbated urban heat stress during nighttime, while during midday both urban and suburban sites experienced comparable heat stress conditions. It is also worth mentioning that the diurnal pattern of UHSI during HWs resembled to some extent the typical UHII pattern with maxima during nighttime and minima during daytime.
According to UTCI, the magnitude of UHSI revealed strong diurnal variation in non-HW conditions, with positive and higher values up to 4.2 • C occurring during late evening and nighttime hours, suggesting in general higher heat stress at the urban site (Figure 7b). UHSI responded to strong HWs by amplification during nighttime, when its magnitude reached 6.4 • C (being higher by at least 2.0 • C compared to non-HW conditions) and decreased during midday (11:00-16:00 LST). The mean diurnal pattern of the UHSI during the strong HW presented a kind of U-shaped curve, though more flattened than the corresponding steeper curve based on HD. This difference is probably attributed to the impact of insolation and wind on thermal stress, since these variables were incorporated in the UTCI calculations but they were omitted in the HD calculations.
It is noticeable that the diurnal pattern of the UHSI difference between HW and non-HW periods was very similar for the two bioclimatic indices (Figure 7a,b), demonstrating a robust intensification of nighttime heat stress conditions in urban, compared to non-urban, sites during severe HWs. Nevertheless, strong HW conditions attenuated the amplitude of UHSI during midday, almost equalizing the heat stress conditions between urban and non-urban sites. This is not surprising given the fact that severe HWs are commonly associated with large scale atmospheric anomalies and persistent anticyclonic conditions that predominate over more local or mesoscale phenomena like UHIs, and this was particularly the case for the extreme HW of 2007 in Greece [34,59]. Although bioclimatic indices involve additional atmospheric variables apart from the air temperature, the observed findings are in qualitative agreement with relevant studies based on the standard metric of UHI intensity. For example, exceptionally hot conditions in Athens have been found to exacerbate significantly the nocturnal UHI intensity (based on the daily minimum air temperature) in contrast to the daytime UHI intensity (based on the daily maximum air temperature) which remained almost constant or even declined with respect to normal summer conditions [59]. Other relevant studies also report stronger synergistic interactions between UHIs and HWs during nighttime, possibly related to increased heat storage and anthropogenic heat release (use of cooling systems) at the urban versus non-urban areas [24,25,27]. Other key contributors for the interaction between UHIs and HWs concern different evaporating rates and unsymmetrical changes in latent and sensible heat fluxes between urban and non-urban sites [24][25][26]. Advective phenomena and local circulations like sea breezes along with increased evaporation rates at coastal non-urban stations have been also found to cause exacerbation of UHI magnitude during HWs, which is more pronounced in daytime hours [26].

Discussion and Conclusions
Heat waves (HWs) represent one of the deadliest natural hazards worldwide, with amplified impacts in urban environments, related to increased population density and the UHI effect. Specific characteristics of HWs such as intensity, duration or seasonal timing combined with population vulnerability have been found to control devastating impacts on human health [2].
The present study focuses on a vulnerable area of the eastern Mediterranean, already characterized as a 'hot spot' with respect to heat-related risk [3,31,32,60] and investigates the change in heat stress levels during HW compared to non-HW (normal summer) conditions as well as the way that heat stress levels respond to HWs in an urban, compared to non-urban, environment.
The adoption of a HW metric accounting for both the intensity and duration of the hot event yielded a total of 46 HW episodes over a nearly 60-year period, but with very rare occurrence until the late 1990s and a profound increased frequency thereafter, reaching up to 5 HWs per year. The duration of the detected HWs ranges between 3 to 11 days while the absolute (AbsExc) and cumulative (CumExc) exceedance of T max reach up to 7.8 • C and 26.7 • C, respectively. These findings for Athens are similar with the results reported in other Mediterranean or eastern Europe cities. For instance, a significant increase in both the number of HW days and the number of long HWs have been recently reported in Madrid, Palma, Nicosia, Rome, Sofia, Zagreb and Ljubljana [19,61,62].
The UTCI and HD indices were used in the study for the assessment of heat stress levels. The CDF curves under HWs and non-HWs conditions reveal a difference of at least one thermal stress category of the index scale based on both indices. In addition, the shift towards higher thermal risk during HWs is apparent across the entire range of the thermal stress distribution. According to HD, the cumulative probability for dangerous conditions that potentially could lead to heat stroke is negligible during non-HW periods, while the corresponding probability during HW periods is 3.6%. Concurrently, the cumulative probability for at least very strong heat stress according to UTCI is extremely low (2.6%) in non-HW periods, while it represents a significant proportion (25.2%) for HW periods. An important issue concerning humans' well-being, especially those living and working in densely urban areas, is the time length in the day when the human body attempts to cool down and recover from heat stress [5,63,64]. The analysis reveals that the cumulative probability for no thermal stress (or slight cold stress) shares a significant proportion (57.3%) of the distribution in non-HW periods, while the percentage is significantly reduced to nearly half (27.3%) during HW periods. The latter implies a considerable shrinkage of the time under no thermal stress during HWs compared to non-HW periods.
Another important finding of the analysis concerns the diurnal variation of thermal stress levels during HWs and non-HWs periods. In particular, the analysis revealed higher diurnal range in the index values during the HWs compared to the non-HW periods. The maximum values of heat stress are observed in midday hours with the value during HWs being 7.4 • C for the UTCI and 7.3 scale points for the HD higher compared to the non-HW periods. Both indices suggest very strong heat stress in the midday hours during HWs, but also discomfort conditions even during nighttime hours.
It is also important to stress, that, even HWs considered to be less severe according to the ranking based on CumExc, may pose strong heat stress to population.
Finally, the study explored the possible interactions or synergies between HWs and heat stress difference between urban and non-urban sites, in accordance to recent research concerning interactions between HWs and UHIs. According to recent literature, there is evidence of synergistic interactions between UHIs and HWs, resulting in the exacerbation of UHI intensity during extreme heat events [24][25][26][27][28]. According to other studies, bigger or smaller cities respond differently during HWs in regards to UHI modification [29]. Other studies have found positive response of heat stress to normal summer conditions between urban and non-urban sites especially during midday hours [65]. The question that arises is, how possible UHI amplification during HWs affects thermal comfort in urban and non-urban sites. In this respect, the indicator of 'Urban Heat Stress Intensity' (UHSI) was applied, defined as the difference between urban and non-urban thermal stress levels as the latter are reflected by the UTCI and HD indices.
Our analysis indicated a strong synergy between HWs and UHSI resulting in marked amplification of UHSI during nighttime hours. The diurnal pattern of the UHSI difference between HW and non-HW periods is very similar for the two bioclimatic indices, demonstrating a robust intensification of nighttime heat stress conditions in urban compared to non-urban sites during severe HWs. Nevertheless, strong HW conditions attenuate the amplitude of UHSI during midday, almost equalizing heat stress conditions between urban and non-urban sites. The magnitude of UHSI based on UTCI (HD) exceeds 6.0 • C (7.0 scale points) during nighttime, being up to at least 2.0 • C (5.0 scale points) larger compared to non-HW conditions.
A crucial question that arises when bioclimatic indices are applied is if these metrics can efficiently capture and imprint the thermal stress in different thermal conditions compared to the actual human thermal responses. Indices have been designed aiming to achieve this, from the simpler to the more comprehensive indices, regardless of whether they achieve it and to what extent. In addition, even the rational indices which take into consideration the mechanisms of heat exchange and the thermo-physiological process on the human body, as well as behavioral factors, do not take into account psychological and expectation factors. The latter factors seem to affect actual thermal responses of people, playing an important role in thermal comfort [66,67]. Nevertheless, recent studies have shown that UTCI satisfactorily predicts the outdoor thermal comfort, based on questionnaires in China [68], in Hong Kong [69] and in Brazil [70]. In the city of Athens, as the area of interest in this research, the surveys of Pantavou et al. [71,72] have shown a strong association between the predicted classes of UTCI and the thermal sensation votes, while UTCI presents the best applicability in terms of predictions and thermal sensation votes among 54 indices. In addition, the thermal sensation votes are predicted fairly well by the UTCI in the heat stress categories of its assessment scale and the UTCI is the index which best simulates the thermal sensation votes among the studied indices [71,72]. Surprisingly, the simple HD index, which is also employed in our study, predicts relatively well the actual human thermal responses and it is ranked first among indices designed for warm weather conditions [71]. Nevertheless, the shortcomings in HD formulation are mainly related to the inability to account for convective and radiative phenomena and the fact that it was initially designed for temperate climates [22]. It is also worthy to add at this point, that adaptation factor is also important in the perception of the thermal environment and thus, the impact of the same weather conditions in a different temperate climate could possibly be more severe [66,[73][74][75].
Our results highlighted the ongoing increasing urban thermal risk, through the increasing frequency of HWs, increasing probabilities for heat stress under HW periods and a positive feedback between HWs and urban heat stress, posing an additional, severe thermal risk to urban population, especially during nighttime.
Funding: This research was funded by the project "THESPIA II-Foundations of synergistic and integrated management methodologies and tools for monitoring and forecasting of environmental issues and pressures" (MIS 5002517) which is implemented under the action "Reinforcement of the Research and Innovation Infrastructure", funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund).